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Vorwort des Herausgebers 


Die Fahrzeugtechnik ist kontinuierlich Veränderungen unterworfen. Klimawan- 
del, die Verknappung einiger für Fahrzeugbau und -betrieb benötigter Rohstoffe, 
globaler Wettbewerb, gesellschaftlicher Wandel und das rapide Wachstum 
großer Städte erfordern neue Mobilitátslósungen, die vielfach eine Neudefi- 
nition des Fahrzeugs erforderlich machen. Die Forderungen nach Steigerung 
der Energieeffizienz, Emissionsreduktion, erhöhter Fahr- und Arbeitssicherheit, 
Benutzerfreundlichkeit und angemessenen Kosten sowie die Möglichkeiten der 
Digitalisierung und Vernetzung finden ihre Antworten nicht aus der singulären 
Verbesserung einzelner technischer Elemente, sondern benötigen Systemver- 
ständnis und eine domänenübergreifende Optimierung der Lösungen. 


Hierzu will die Karlsruher Schriftenreihe für Fahrzeugsystemtechnik einen Bei- 
trag leisten. Für die Fahrzeuggattungen Pkw, Nfz, Mobile Arbeitsmaschinen 
und Bahnfahrzeuge werden Forschungsarbeiten vorgestellt, die Fahrzeugsys- 
temtechnik auf vier Ebenen beleuchten: das Fahrzeug als komplexes, digi- 
talisiertes mechatronisches System, die Mensch-Fahrzeug-Interaktion, das 
Fahrzeug in Verkehr und Infrastruktur sowie das Fahrzeug in Gesellschaft 


und Umwelt. 


Eine automatisierte Fahrzeuglängsführung entlastet den Fahrer und erlaubt da- 
rüber hinaus eine Steigerung der Energie- und Zeiteffizienz, ohne das Situati- 
onsbewusstsein negativ zu beeinflussen. Häufig wird für die vorausschauende 
Optimierung der Geschwindigkeitstrajektorie die Methode der dynamischen 
Programmierung eingesetzt, die das globale Optimum liefert, jedoch sehr hohe 
Anforderungen an Rechenleistung und Speicherplatz stellt. 


Vorwort des Herausgebers 


Herr Jauch schlägt in seiner Arbeit ein rekursives Verfahren zur Generierung 
optimierter Geschwindigkeitstrajektorien unter Nebenbedingungen für beliebig 
lange Strecken vor, dessen Rechenaufwand klein ist und lediglich linear mit der 
Streckenlänge steigt. Es basiert auf der Anpassung einer B-Spline-Funktion 
an kartenbasierte Datenpunkte und ist sowohl für lineare als auch für nichtlin- 
eare Probleme geeignet. Es liefert in der Anwendung auf batterieelektrische 
Pkw eine deutliche Erhöhung der Energieeffizienz bzw. der Durchschnitts- 
geschwindigkeit. 


Karlsruhe, im Juli 2023 Frank Gauterin 


Abstract 


This work describes a novel method for approximating an unbounded number 
of data points using a B-spline function in the linear and nonlinear weighted 
least squares sense. The developed method is based on iterative algorithms 
for state estimation and its computational effort increases linearly with the 
number of data points. The method allows to shift the bounded definition 
range of a B-spline function during run-time such that the latest data point can 
be considered for approximation regardless of the initially chosen definition 
range. Furthermore, the shift operation allows to decrease the sizes of matrices 
in the state estimators in order to reduce computational effort in both offline 
applications, in which all data points are available at once for processing, and 
online applications, in which additional data points are observed in each time 
step. 


The trajectory optimization problem is formulated such that the approximation 
method computes a B-spline function that represents the desired velocity trajec- 
tory with respect to time using data points created from map data. The compu- 
tational effort of the resulting direct trajectory optimization method increases 
only linearly with the unbounded temporal length of the planned trajectory. The 
combination with an adaptive model that describes the power train properties 
of a battery electric vehicle with fixed gear box transmission ratio allows to 
optimize velocity trajectories with respect to travel time, comfort and energy 
consumption. 


The trajectory optimization method is extended to a driver assistance system for 
automated vehicle longitudinal control that is tested in simulations as well as in 
real test drives. Simulated drives on a chosen reference route need up to 3.4 % 
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Abstract 


less energy with the automated longitudinal control than with a human driver 
at the same average velocity. For the same energy consumption the automated 
longitudinal control achieves a 2.6 % higher average velocity than a human 
driver. 
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Kurzfassung 


Diese Arbeit beschreibt ein neuartiges Verfahren zur linearen und nichtlinearen 
gewichteten Kleinste-Quadrate-Approximation einer unbeschränkten Anzahl 
von Datenpunkten mit einer B-Spline-Funktion. Das entwickelte Verfahren 
basiert auf iterativen Algorithmen zur Zustandsschätzung und sein Rechen- 
aufwand nimmt linear mit der Anzahl der Datenpunkte zu. Das Verfahren 
ermöglicht eine Verschiebung des beschränkten Definitionsbereichs einer B- 
Spline-Funktion zur Laufzeit, sodass der aktuell betrachtete Datenpunkt un- 
geachtet des anfangs gewählten Definitionsbereichs bei der Approximation 
berücksichtigt werden kann. Zudem ermöglicht die Verschiebeoperation die 
Reduktion der Größen der Matrizen in den Zustandsschätzern zur Senkung des 
Rechenaufwands sowohl in Offline-Anwendungen, in denen alle Datenpunkte 
gleichzeitig zur Verarbeitung vorliegen, als auch in Online-Anwendungen, in 
denen in jedem Zeitschritt weitere Datenpunkte beobachtet werden. 


Das Trajektorienoptimierungsproblem wird so formuliert, dass das Approxi- 
mationsverfahren mit Datenpunkten aus Kartendaten eine B-Spline-Funktion 
berechnet, die die gewünschte Geschwindigkeitstrajektorie bezüglich der Zeit 
repräsentiert. Der Rechenaufwand des resultierenden direkten Trajektorienop- 
timierungsverfahrens steigt lediglich linear mit der unbeschränkten zeitlichen 
Trajektorienlänge an. Die Kombination mit einem adaptiven Modell des 
Antriebsstrangs eines batterie-elektrischen Fahrzeugs mit festem Getriebeüber- 
setzungsverhältnis ermöglicht die Optimierung von Geschwindigkeitstrajekto- 
rien hinsichtlich Fahrzeit, Komfort und Energieverbrauch. 


Das Trajektorienoptimierungsverfahren wird zu einem Fahrerassistenzsystem 
für die automatisierte Fahrzeuglängsführung erweitert, das simulativ und in 


Kurzfassung 


realen Erprobungsfahrten getestet wird. Simulierte Fahrten auf der gewähl- 
ten Referenzstrecke benötigten bis zu 3,4 % weniger Energie mit der auto- 
matisierten Längsführung als mit einem menschlichen Fahrer bei derselben 
Durchschnittsgeschwindigkeit. Für denselben Energieverbrauch erzielt die au- 
tomatisierte Längsführung eine 2,6 % höhere Durchschnittsgeschwindigkeit 
als ein menschlicher Fahrer. 
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1 Introduction 


1.1 Background 


The topic of this work falls into the context of driver assistance systems for 
automated vehicle longitudinal control based on map data. Such systems can 
increase safety, improve comfort and reduce the energy consumption of the 
vehicle by computing an energy-efficient course of velocity for the road section 
ahead. 


In the recent past, energy savings have mainly been achieved by power train 
improvements such as downsizing of the internal combustion engine as well 
as hybrid and purely electric power trains. When taking into account the in- 
creasing costs for further power train improvements, assistance systems for 
energy-efficient longitudinal control offer a more cost-efficient way to obtain 
additional energy savings. According to [117, p. 141], efficient driving has the 
second largest potential for energy savings (30 %) after electrification of the 
power train (75 %). 


Driver assistance systems for automated longitudinal control (ALC) compute 
the desired course of velocity, also called velocity trajectory, by solving a 
trajectory optimization problem. Different optimization approaches are known, 
e.g. Dynamic Programming (DP) and direct methods (DM). 


The popular DP approach finds the global optimum for a given optimization 
problem and its computational effort grows linearly with the temporal length of 
the planned trajectory. However, only coarsely discretized problems with few 
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state dimensions can be solved in real-time on a series electronic control unit 
(ECU) with limited resources. 


A possible way is to compute a rough long-term trajectory with DP that serves 
as a reference trajectory for a DM. In addition to discretizing the optimization 
problem, DM approximate the system state and control between the discretiza- 
tion steps by a function with respect to time. DM are, however, only suitable 
for short-term trajectories since the effort of DM increases exponentially with 
the temporal trajectory length. 


This work addresses the need for computationally less demanding trajectory 
optimization algorithms for ALC of a battery electric vehicle (BEV). 


1.2 Approach 


The research problem is approached by developing a novel method for the 
generic task of approximating an unbounded number of data points using a B- 
spline function in the linear and nonlinear weighted least squares (WLS) sense. 
The developed method is based on iterative algorithms for state estimation. 


The approximation problem is reformulated as a trajectory optimization prob- 
lem such that the approximation method computes a velocity trajectory with 
respect to time using data points created from map data. The novel trajectory 
optimization method falls into the category of DM and its effort increases only 
linearly with the trajectory length instead of exponentially. The combination 
with an adaptive model that describes the power train properties of the BEV 
allows to plan velocity trajectories whose resulting energy consumption varies 
depending on the chosen relative weighting of different target criteria. 


The trajectory optimization is extented to a driver assistance system for ALC 
that is tested in simulation as well as in real test drives. In simulations on a 
chosen reference route the ALC is compared to a recorded and re-simulated real 
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Figure 1.1: Structure of this work illustrated by the V model that is used to describe the software 
development process. 


drive with manual longitudinal control (MLC) regarding energy consumption 
and average velocity. 


1.3 Outline 


Figure 1.1 illustrates the structure of this work, which follows the V model. The 
V model describes the software development process and distinguishes different 
process levels. On the descending branch of the V model, the development goal 
is analysed and iteratively broken down into subgoals. Thereby requirements 
are specified for subsystems that ultimately need to be implemented. 
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Chapter 2 of the work corresponds to the descending branch of the V model 
and will provide literature background as well as identify research gaps. The 
development goal on system level is a driver assistance system for energy- 
efficient ALC of a BEV (Section 2.1) which is broken down into a trajectory 
optimization method on component level (Section 2.2) and B-spline approxima- 
tion methods on algorithm level (Section 2.3). Section 2.4 gives an introductive 
overview of adaptive filters for Chapter 3 and Chapter 5 without identifying a 
research gap. 


The ascending branch of the V model includes testing and iterative integration 
of subcomponents into larger systems until the main goal is reached. Chapter 3 
to Chapter 6 address this phase by presenting contributions for closing the 
research gaps on each level of the V model. 


Chapter 3 presents novel algorithms for iterative linear and nonlinear WLS 
approximation of data with a B-spline function. The algorithms are referred to 
as recursive B-spline approximation (RBA) and nonlinear recursive B-spline 
approximation (NRBA), respectively. 


Chapter 4 provides neccessary foundations for Chapter 5 and Chapter 6. With- 
out having a direct counterpart to the left branch, it states models of the vehicle 
that are used for determining energy-efficient trajectories and for evaluating 
the energy consumption of the research vehicle in simulations on a chosen 
reference route. 


Chapter 5 applies RBA and NRBA to the trajectory optimization problem and 
presents extensions for taking into account trajectory constraints. 


Chapter 6 describes the development process of the ALC and its system archi- 
tecture including the interaction of its components. Furthermore, Chapter 6 
investigates the energy-saving potential with the ALC on the reference route 
in comparison to a recorded and re-simulated real drive with MLC as well as 
effects of ALC parameters on the energy consumption that results with the 
ALC. 


1.4 Works, contributions and support by others 


Chapter 7 concludes the journey through the V model, summarizes the contri- 
butions to research gaps and gives suggestions for future research. 


1.4 Works, contributions and support by 
others 


The research described in this work took place within the research project 
e-volution [41] funded by the German Federal Ministry of Education and Re- 
search (support code 16EMO0071). The Dr. Ing. h.c. F. Porsche AG was 
among the research partners and supported the project work by providing a 
battery electric research vehicle for test drives as well as vehicle data and mea- 
surement data. Furthermore, an employee of the Porsche Engineering Services 
GmbH integrated the developed ALC system into the prototyping ECU of the 
research vehicle and conducted test drives. 


The provided research vehicle was developed during the preceding research 
project e-generation [22], which was also funded by the German Federal Min- 
istry of Education and Research (support code 16N 11865) and realized in 
cooperation with the above mentioned companies. The author of this work did 
not contribute to e-generation but reused the following results: 


Despite differing strongly in detail, the B-spline data approximation algorithms 
stated in Chapter 3 adopt the iterative, filter-based characteristics of the previous 
method for polynomial data approximation [22, pp. 29-35], which was patented 
by F. Bleimund and S. Rhode [21]. A detailed differentiation between the 
approaches is given at the end of Chapter 3. 


The adaptive traction force model (ATFM) stated in Subsection 4.7.1 stems 
from F. Bleimund [22, pp. 24-29] and is applied without modifications. The 
idea to model the power train characteristics with a kernel regression model, as 
done in Subsection 4.7.2 using the adaptive electrical power model (AEPM), 
originates from S. Rhode [146], who also provided a first script in MATLAB 
that performs this task. Based on this, both authors continued the investigations 
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independently from each other. Problem-specific adaptions for using the model 
for trajectory optimization as well as the determination of model hyperparame- 
ters were done by the author of this work. 


Trajectory optimization within the previous ALC [22, pp. 29-35] developed by 
F. Bleimund differs significantly from the one described in this work, although 
both are examples of DM. The previous ALC defined trajectories with respect 
to position by a cubic polynomial. The few degrees of freedom of this function 
only allowed to represent trajectories with simple shape over short distances. 
The energy that the vehicle will need for tracking the planned trajectories was 
not considered in the optimization process. 


Therefore, regarding trajectory optimization this work only adopts the genera- 
tion of an upper speed limit (Section 5.1) from e-generation [22, pp. 29, 30]. 
The adaptions made are limited to parameterization improvements. At an early 
stage, with non-final trajectory representation, AEPM and nonlinear filter, tests 
for finding suitable target criteria weightings for the trajectory optimization 
including the elecrical power (Section 5.4) were done by A. Thorgeirsson [173] 
as a student assistant. 


Moreover, reused results include agreed interfaces between the ALC system 
and the vehicle, the technical architecture of the ALC system as well as a 
Hardware-in-the-Loop (HiL) test bench, which simulates the research vehicle 
on the reference route and runs the developed ALC on the same ECU type as 
in the real research vehicle. The HiL test bench (Section 6.1) was created by B. 
Fath [22, pp. 60-73]. 


The route data module (Subsection 6.2.1) was implemented by D. Dorr [22, 
pp. 23, 24]. Within this work, the parameterization was enhanced and logic 
for processing roundabout data was added. The parameter adaption module 
(Subsection 6.2.2) including ATFM was developed by F. Bleimund [22, pp. 24- 
29]. The author of this work added the AEPM. 


1.5 Prepublications and their citations 


To the controller module the author of this work added the functionality stated 
at the beginning of Subsection 6.2.4 and replaced the original proportional- 
integral-derivative (PID) control by a model predictive control (MPC). Overall 
architecture and control loop including the pilot control (Subsubsection 6.2.4) 
were adopted from F. Bleimund [22, pp. 35, 36]. First implementations and 
simplified tests of the MPC (Subsubsection 6.2.4) were done by C. Lee [98] 
during his time as student assistant after completing his Master’s thesis. 


First evaluations of the energy-saving potential (Section 6.3) with non-final 
ALC setup were conducted by A. Thorgeirsson [173] while being employed as 
a student assistant. For more detailed differentiations regarding the individual 
contributions, the reader is referred to the afore-stated parts of this work. 


1.5 Prepublications and their citations 


The research community benefits from fast publication of new findings. There- 
fore, parts of this dissertation have been prepublished. Namely, RBA and 
NRBA described in Chapter 3 are the topics of [78, 80] and the contents of 
Chapter 4, Chapter 5 and Chapter 6 are mentioned in a very condensed way in 
[78, 79]. 


This work will put the topics of these publications into a broader context and 
present them in much more detail. For more flexibility in rearranging content 
compared to a cumulative dissertation, only parts of previous publications have 
been added. To avoid interruptions of the reading flow because of frequent 
changes between paraphrases and direct quotes, all content adopted from pre- 
publications, regardless of the extent of adaptions made, i.e. both paraphrases 
as well as direct quations, will be indicated as follows: 


o ... Prepublished content... Adopted from [prepublication]. 
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For reproducibility and convenient further use by other researchers the MATLAB 
source code for RBA, NRBA and trajectory planning has been provided along 
with the corresponding publication [75, 76, 77]. However, the source code of 
the trajectory planning method was published with a generic vehicle model 
instead of the parameters of the research vehicle and for the simplest planning 
case only. Therefore the results presented in Chapter 4 to Chapter 6 are not 
exactly reproducible by others. 


1.6 Notation 


In this work vectors are printed in bold font and matrices are printed bold in 
capital letters. For example, a and A are scalars, a is a vector and A is a matrix. 
Unless explicitly stated, the sizes of vectors and matrices are assumed to be 
chosen appropriately. 


A hat above a variable a such as â indicates an estimate of a or the solution of 
an optimization problem with respect to the variable a. 


This work also adopts the MATLAB-like index notation. The first index of a 
matrix A refers to rows, the second to columns. For example, A; ; refers to the 
entry in row i and column j of matrix A. The colon operator (:) refers to all 
rows or columns, or a range. A. ; means all rows of column j. Analogously, A; 
denotes all columns of row i. A; is an abbreviation of A;,.. The colon operator 
in the index of Aj.4 ; extracts a column vector of rows one to four in column j 
of matrix A. Additionally, a[i] refers to the i-th element of vector a. 


| - | denotes the absolute value and ()" the transpose operation. 


2 Scientific and technical state of 
the art 


This chapter provides background information and identifies research gaps 
regarding driver assistance systems for ALC in Section 2.1, trajectory opti- 
mization in Section 2.2 and types of spline representations as well as B-spline 
approximation methods in Section 2.3. Furthermore, Section 2.4 gives an 
overview on adaptive filters for Chapter 3 and Section 4.7. 


2.1 Assistance systems for automated 
longitudinal control 


The tasks that the driver performs can be divided into three levels: 


e Navigation: Choosing a travel route 
e Guidance: Determining a suitable speed from the traffic environment 
e Stabilization: Aligning the vehicle motion with the guidance variables 


Driver assistance systems can support on each of these levels. 


O Regarding the vehicle guidance, assistance functions can be distinguished by 
their operating mode as depicted in Table 2.1 [187, pp. 36-40]: 


e Informing and warning functions 
e Continuously automated functions 
e Intervening emergency functions 
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Table 2.1: Operating modes regarding vehicle guidance [187, p. 38, adapted] 


Mode Informing and warning Continuously automated Intervening emergency 
name functions functions functions 
Mode A B C 

character 


Vehicle Only indirect control via + Immediate direct control Immediate direct control 


control the driver « Division of driving tasks in near-accident 
between human driver and situations or situations 
function that cannot handled by an 
e Control always remains average driver 
overridable 
Examples + Display of current speed Usually convenience Usually safety functions: 
limit based on traffic sign functions: Automatic emergency 
recognition + Adaptive cruise control braking (system 
e Lane departure warning + Lane centering triggered) 


This work deals with continuously automated functions for longitudinal control. 
Such functions can contribute to increasing safety and comfort as well as to 
reducing energy consumption to different extents. 


For example, the cruise control (CC) function controls the brake and engine 
torque so that the vehicle maintains the speed selected by the driver. Adaptive 
cruise control (ACC) is an extension of CC that uses a radar sensor to detect 
a vehicle ahead. ACC can reduce the selected vehicle velocity to maintain a 
chosen time gap to a vehicle ahead [52, 196]. Studies cited in [187, pp. 1140- 
1145] report that with activated ACC drivers mention feeling safer and more 
relaxed compared to manual driving. Adopted from [79]. ACC-like systems 
increase comfort by reducing the driver workload. However, the effects on 
safety are not as clear as the effects on comfort and require more comprehensive 
investigations including psychological effects. These are adressed by works that 
are mentioned in the above cited literature. Briefly summarized, automation of 
tasks favors that the driver diverts the attention to things that are irrelevant for 
driving. In assistance systems that still require monitoring by the driver, this 
can cause a Safety issue if the driver unexpectedly needs to take back control. 
For example, in one of the studies, drivers looked away from the traffic ahead 
for longer periods when driving with ACC [187, pp. 1140-1145]. 
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2.1 Assistance systems for automated longitudinal control 


Table 2.2: Time and energy savings with automated longitudinal controls that consider map data 
compared to benchmark control methods. Dynamic Programming (DP), direct methods 
(DM), conventional vehicle (CV), hybrid electric vehicle (HEV), battery electric vehicle 
(BEV), urban road (U), country road (C), highway (H), artificial road profile (A). 


Source Optimization Vehicle Energy difference Time difference Benchmark Road 
[114 DP CV -10 % 0% Human U,C 
[139 DP CV -23.8 % 1.5 % Human U,C 
[139 DP CV -10.2 % -1.3 % Human U,C 
[100 DP CV -12 % -3 % Human UG 
[100 DP CV -8.28 % 0.01 % Human U,C,H 
[100 DP CV 1.71% -3.60 Yo Human U,C,H 
[183 DP HEV -18.1 % -1.15 % Heuristic strategy U,C,H 
[183 DP HEV 0 % -21.13 % Heuristic strategy U,C,H 

[86] DM HEV —4 % n/a Constant velocity H 
[188 DP BEV —5.83 % n/a PID control A 
[188 DM BEV 5.40 % n/a PID control A 
[188 DM BEV -1.27 % n/a PID control H 
[188 DM BEV -1.14 % n/a PID control H 


Driving efficiency oriented extensions of ACC use a corridor of allowed dis- 
tances to the vehicle ahead instead of a specific distance value given by the 
selected time gap. Within this corridor the vehicle is controlled according to 
an optimization that considers the energy consumption of the vehicle, either 
explicitly with an energy consumption model [18, 88, 127, 195] or implicitely 
using a proxy for the energy consumption, such as acceleration [1, 107, 112]. 
The effectiveness of both approaches is compared in [81]. Section 4.7 will 
address consumption models in more detail. 


O One step further go ACC enhancements that additionally determine the ap- 


propriate speed depending on map data so that the vehicle automatically slows 
down if a curve is ahead. Map information enables such systems to determine 
an even more energy-efficient or time-efficient driving strategy as Table 2.2 
illustrates. Adopted from [79]. 


Further developments of assistance systems increase the degree of automation 
and finally lead to autonomous vehicles that do not require a driver. The degree 
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Table 2.3: SAE J3016 levels of driving automation [172, adapted]. 


Level name No Driving Driver Partial Driving Conditional High Driving Full Driving 
Automation Assistance Automation Driving Automation Automation 
Automation 
Level number 0 1 2 3 4 S 
Type Driver support system Automated driving system 
Vehicle Human is driving when the system is Human is not driving when the system is 
control engaged, even if pedals and steering wheel engaged, even if sitting in the driver’s seat 
are not touched 
Monitoring Human must constantly supervise the system On request by Feature will not require 
and fallback and override it as the situation requires the system, human to take over driving 
solution human must 
take over 
driving 
Capability Provide Provide Provide Can drive the vehicle under Can drive the 
warnings and EITHER longitudinal limited conditions as longas vehicle under 
momentary longitudinal AND all required conditions are met all conditions 
assistance OR lateral control 
lateral control support 
support 
Example Display of EITHER Lane centering Traffic jam Local Same as 4 but 
current speed lane centering AND chauffeur driverless taxi system can 
limit based on OR adaptive cruise drive 
traffic sign adaptive cruise control Installation everywhere 
recognition control SIMULTANE- of pedals and in all 
OUSLY steering wheel conditions 
Lane optional 
departure 
warning 
Automatic 
emergency 
braking 


of automation according to SAE International is depicted in Table 2.3. Both 
mode A and mode C functions according to Table 2.1 are SAE level 0 functions 
because they are only capable of warnings and momentary assistance. Most 
ACC systems fall into the category SAE level 1, because they can perform 
longitudinal control in certain use cases but the driver has to monitor them, 
has to override them occasionally and needs to perform the lateral control 
by operating the steering wheel. Some system enhancements include lateral 
guidance [187, p. 1146]. Systems that offer both longitudinal and lateral control 
belong to at least SAE level 2. 
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O ALC systems that determine an energy-efficient driving strategy mainly have 
been developed for vehicles with internal combustion engine [139] or hybrid 
electric vehicles [183]. The degrees of freedom of the power train in a con- 
ventional vehicle include motor torque, clutch state and selected gear. The 
various operating modes of a hybrid power train translate to additional degrees 
of freedom. Therefore an energy-efficient driving strategy is the solution to a 
nonlinear and high dimensional problem. To perform the required computa- 
tions on a ECU with limited computational power in real-time, the referenced 
approaches need to reduce the problem complexity, discretize it coarsely and 
partly use parallelizable, iterative and approximative approaches. 


Research gap: In contrast to conventional and hybrid vehicles, aBEV usually 
has a 1-speed gear box and therefore a constant gear ratio. If the BEV has 
multiple motors, as the research vehicle does, they have at least similar char- 
acteristics so that the possible benefits from including the torque distribution 
in the optimization problem over a rule-based torque distribution strategy that 
considers efficiency and driving stability are small. With a rule-based strategy, 
as described in Section 4.3 for the research vehicle, the only degree of freedom 
that is left in the power train for optimization is the total torque. The power 
required by auxilaries can also be significant but will not be considered in the 
optimization based on the rationales mentioned in Section 4.4. 


The lack of ALC systems that take advantage of this resulting simple structure 
by incorporating a trajectory optimization approach that is on the one hand 
less suitable for the complexity of the power train of a conventional or hybrid 
vehicle but on the other hand oriented towards low computation power demand 
by design, such as a local optimization based on adaptive filters, an overview of 
which is given in Section 2.4, poses a research gap. The next section will illus- 
trate and compare the features of different trajectory optimization approaches 
in more detail. Adopted from [78, 79]. 
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2.2 Trajectory optimization 


Driver assistance systems that perform continuously automated functions cal- 
culate input commands to steering, brakes, engine and power train in order to 
achieve a desired vehicle motion, also called trajectory. Trajectory planning is 
also known as motion planning in robotics [187, p. 1414] and differs from a 
path planning in that it additionally assigns a time law to a geometric path [55]. 
Reviews of state of the art motion planning methods for automated driving are 
provided in [40, 59, 90]. 


With increasing extent of automation and number of degrees of freedom the 
decision making process of driver assistance systems becomes more complex 
and the trajectory is found by solving an optimization problem. Trajectory 
optimization denotes the process of determining a trajectory of a dynamical 
system including the control input to the system. Thereby the trajectory must 
meet the system constraints and optimize a performance measure [185, p. 3]. 
This section gives an overview of trajectory optimization methods. 


The trajectory is often optimized with respect to a trade-off between target 
criteria such as comfort, safety, energy comsumption and travel time. Trajectory 
constraints result from the vehicle dynamics, for example power restrictions, 
as well as from the environment, for example lanes and other vehicles [187, 
p. 1414]. 


In static optimization problems the optimal values of a finite set of variables p 
need to be determined such that a cost function J(p) is minimized. Trajectory 
optimization problems fall into the category of dynamic or infinite dimensional 
optimization because the optimization variables are functions x(t) of an inde- 
pendent variable r, usually time. Assessing the performance of the trajectory by 
a scalar quantity requires a cost functional, which is a function of a function. 
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Optimal control problem 


Dynamic 
Programming 


Indirect methods Direct methods 


Figure 2.1: Methods for solving the optimal control problem 


The trajectory optimization problem can be mathematically formulated as an 
optimal control problem (OCP): 


u*(t) = argmin J(u(t)) 
u(t) 
J(u(t)) = ie (x(t), u(t), t) dt + V (xte), tf) 
(2.1) 
subject to 


X(t) = f (x(t), u(t),t), x(0) = xo, 
g(x(tp).tp) =0, h(x(0),u(0),1) < 0 Vr € [0,17]. 


In an OCP we seek for t € [0, tf] the control input trajectory u(t) € R” for 
a system with state x € R” that leads the system from its initial state xo to 
a terminal state x,, while minimizing the cost functional J and meeting the 
system dynamics model f, equality constraints g and inequality constraints 
h. u*(t) is the optimal input trajectory and x; the resulting state trajectory. J 
includes integral costs / and terminal costs V [187, pp. 1415-1416]. 


O Figure 2.1 illustrates that known solution approaches to the trajectory opti- 


mization problem can be categorized into the three principles Dynamic Pro- 
gramming (DP), indirect methods (IM) and direct methods (DM). Adopted 
from [78, 79]. However, some approaches combine characteristics of several 
principles. Therefore the assignment is not necessarily unique [131, p. 27]. 
Apart from simple problems and special cases, solutions are derived numeri- 
cally [185, p. 9]. The following subsections describe the three principles along 
with application examples. 
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Dynamic Programming 


For purely continuous systems a solution to the Hamilton-Jacobi-Bellman par- 
tial differential equation provides sufficient optimality conditions for an optimal 
solution of the control problem. The DP algorithm samples the continuous state 
and control spaces in order to achieve a discrete OCP: 


u*(t) = argmin J(u(t)), J(u(t)) = > I(x(k),u(k),k) 
u(t) k=0 


(2.2) 
subject to 


x(k +1) =f (x(k), u(k)), x(0) = xo, k =0,...,k¢ -1 


The optimal control sequence u*(k) that minimizes the sum of cost from the 
initial state to the final state is computed iteratively using value iteration. Equal- 
ity or inequality constraints are taken into account by setting the costs / = oo 
for cases that violate constraints [187, pp. 1425-1430]. 


The algorithm takes advantage of the optimality principle of Bellman, which 
states that an optimal trajectory is composed of optimal subtrajectories. This 
principle allows to split the OCP into many smaller problems. 


Since the value function depends on the system state x, the optimal control 
is a closed feedback control loop. The closed-loop control is valid within the 
bounds of the sampled state space [131, p. 5]. 


DP is especially beneficial for trajectory planning in driving scenarios with 
dynamic environment and constraints resulting from other traffic participants 
[87] and for planning parking maneuvers [103]. DI An ACC based on DP 
is proposed in [53]. DP based algorithms for energy-efficient ALC exist for 
vehicles with internal combustion engine [100, 139], hybrid electric vehicles 
[183], plug-in hybrid electric vehicles [197] and BEVs [188]. Adopted from 
[78, 79]. 
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If the time when a vehicle reaches a location is unimportant, the problem can 
be simplified. However, without temporal dependency the decision graph is 
cyclic. Therefore the value iteration needs to be replaced with another graph 
searching method that draws on the DP paradigm, for example the Dijkstra 
algorithm or the faster A* algorithm [187, p. 1430]. In the DARPA urban 
challenge 2007 many successful autonomous vehicles applied an A* based path 
planning algorithm [34]. 


Indirect methods 


Indirect methods (IM) transform the OCP into the problem of solving a system 
of nonlinear equations. First, they determine the neccessary first-order optimal- 
ity conditions for an OCP, which leads to a boundary value problem consisting 
of a set of differential equations, the so-called Hamilton equations. Then they 
solve the boundary value problem numerically, e.g by a Newton method [185, 
pp. 10, 11]. 


For example, variational calculus requires that the first derivative of the func- 
tional J(u(r)) diminishes for the optimal control u*(t). The Lagrange mul- 
tiplier method allows incorporating differential equality constraints resulting 
from the system dynamics. Without inequality constraints the Hamilton equa- 
tions then read 


N . ol faf] 2 ôl A 
x = f(x,u,t), O [| A, A LA A. (2.3) 


Herein A denotes the Lagrange multipliers [187, pp. 1417-1418]. The approach 
1s also known as maximum principle of Pontryagin [131, p. 32] or minimum 
principle of Pontryagin [185, p. 10], depending on the formulation of the opti- 
mization problem. 


IM only state a necessary condition for optimality and considering state con- 
straints is often difficult [127]. Furthermore, they require initial conditions 
for the Lagrange multipliers. This made them unsuitable for many automotive 
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applications [187, p. 1420]. Nevertheless, there are indirect trajectory optimiza- 
tion applications in literature, e.g. for conventional vehicles in [126, 160] and 
for electric vehicles in [39, 132]. 


Direct methods 


o DM compute the solution of an OCP for a continuous system by discretizing 
it with respect to time. Between the discrete time steps they approximate the 
system state and control with a function of time. Thereby DM translate the 
infinite dimensional problem into a finite dimensional optimization problem. 
Adopted from [79]. System dynamics enter the optimization problem via 
constraints for discrete points in time. Therefore constraints are also finite. 


For example, the input trajectory can be represented by a polynomial ı and its 
finite dimensional parameter vector u: 


u(t) = W(t, ü) (2.4) 
Then the system dynamics are 
X(t) = f(x), y (1,4), 1), x(to) = xo. (2.5) 


This initial value problem can be solved by a differential equation solver. The 
system trajectory is denoted by 


x(t) = p(1,u). (2.6) 
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The resulting static problem that approximates the dynamic OCP reads 
u*(t) = argmin J(ü), 
ü 
tf 
J(ü) = T 1 ($t, @), W(t, ā),t) dt + V (6(,, @)) on 


subject to 


g(dlt,ü),ty) =0, h(dlt,u),y(t,ü))<0, i=1,...,N. 


O The cost functional J is optimized directly using an optimization method 


that varies the finite state and control values of the functional representation. 
Usually the optimization problem is nonlinear and solved using sequential 
quadratic programming methods (e.g. [65, 92]) or interior point methods. 


DM often occur in literature as MPC, which solves (2.7) on a receding horizon 
[187, pp. 1420-1425,1431-1432]. Adopted from [78]. For example, [56] 
applies MPC for generating locally optimal trajectories, [28] proposes a MPC 
based lateral control that tracks the lane centerline while avoiding collisions, 
[47] calculates the control inputs with MPC so that the vehicle tracks a planned 
path and [127] uses MPC for an ACC system that additionally aims to reduce 
the energy consumption. 


Comparison 


A comparison of the stated solution principles is presented in [131, p. 8] and 
further summarized in Table 2.4. The generalized statements apply to many 
methods found in literature but individual methods might differ. 


DP computes a globally optimal control because of the sufficient optimality 
conditions of the Hamilton-Jacobi-Bellman equations. The computational effort 
grows linearly with increasing time horizon and number of discrete states. DP 
offers a closed-loop feedback control that is beneficial for the system stability. 
DP algorithms converge globally in the discrete state space. However, the 
required discretization leads to an exponential increase of computational effort 
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Table 2.4: Comparison of methods for solving an optimal control problem (OCP). The table shows 
average values for common algorithms applied to the same OCP. Some methods may 
not fulfill all of the generalized statements. Summarized version derived from [131, 
p- 8]. *: According to [131, p. 8], an exception that achieves polynomial increase using 
a convexification approach is stated in [151]. 


Dynamic Indirect Direct 

Category i 

Programming (DP) methods (IM) methods (DM) 
Optimality of solution Global Local Local 
Domain of convergence Global Smaller than DM Larger than IM 
Ease of initialization Higher than DM Lower than DM Medium 
Classes of solvable systems Less than IM Medium More than IM 
Control Closed-loop Open-loop Open-loop 


Increase of complexity 


with increased 


- Accuracy Polynomial Linear Polynomial 
- Time horizon Linear Polynomial Exponential* 
- # of continuous states Exponential Polynomial Polynomial 
- # of discrete states Linear Polynomial Polynomial 
Required system knowledge High Higher than DP Lower than DP 


with the number of continuous states. Therefore, DP usually can only be applied 
to systems with a state and control dimension up to four and requires coarse 
discretization. Therefore the solution accuracy is comparatively low. If the 
discretization is refined for more accurate solutions, the computational effort 
increases polynomially. DP requires a high system knowledge and can solve 
less classes of OCP than other principles. 


IM are very accurate and the computational effort only grows linearly when the 
accuracy is increased. With respect to the time horizon and number of states, 
the complexity grows polynomially. IM can solve more OCP classes than DP 
but less than DM. IM only compute locally optimal open-loop controls which 
means that the control loop is not closed using system feedback information. 
IM converge within a smaller domain than DM and in general, applying IM 
requires more knowledge than other methods. 
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Table 2.5: Suitability of Dynamic Programming (DP), direct methods (DM), indirect methods (IM) 
and combinations for selected problem features. Particularly suitable (++), suitable (+), 
less suitable (—). Derived from [187, p. 1430]. 


Many Continuous Global Long 


Approach . . 
states states optimum horizon 
DP ad = + + 
DM + + = = 
DP + DM + + $ + 
DP+DM+IM + + + ++ 


o DM converge within a larger domain than IM. They can deal with numerous 
states, require less control knowledge and can solve more OCP classes than 
other methods. Like IM, DM provide only locally optimal open-loop control. 
The acceptable solution accuracy mainly results from discretization errors. The 
computational effort grows polynomially with the solution accuracy and the 
number of states. With respect to the time horizon, effort usually grows expo- 
nentially because in each time step all controls have to be analyzed. Therefore 
the optimization horizon is mostly restricted to few seconds [131, pp. 2-12, 27- 
37], [187, pp. 1415-1431]. 


Due to their complementary properties, different approaches are combined for 
solving difficult, farsighted trajectory optimization problems. Then DP provides 
a rough long-term plan or reference trajectory for a DM that computes feasible 
trajectories within a short time horizon. Adopted from [78, 79]. An example 
is the energy-optimal adaptive cruise control in [186]. Closed-form solutions 
of calculus of variations can be applied to speed up DP. Table 2.5 shows the 
benefits of combining different methods [187, pp. 1430-1431]. 


Research gap: Both DP and DM are popular for automotive applications on 
their own and combining their orthogonal features offers great potential. How- 
ever, the exponential growth of computational effort with increasing time hori- 
zon limits the application of DM to short time horizons. Hence, the research 
gap regarding trajectory optimization consists of available DM with lower com- 
plexity. 
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2.3 B-spline data approximation 


Data approximation can be viewed as a special case of trajectory optimization 
using the DM approach and instead of a polynomial, a B-spline function can 
represent the input trajectory from (2.4) or the system trajectory from (2.6). 
Following an overview of spline representations, this section states methods 
for determining a data approximating B-spline function and identifies the cor- 
responding research gap. 


Spline representations and their features 


A function that is used for approximating data points needs to possess a suf- 
ficient number of parameters which translates to the number of degrees of 
freedom. A complicated relationship in the data usually requires more function 
parameters in order to achieve an acceptable representation. 


With polynomials the degree is coupled to the number of parameters and as 
the degree is increased, the polynomials become computationally expensive to 
evaluate [5]. Assume a polynomial p(u) given by 


d 
p(u) = Br uc. (2.8) 
i=0 
With a high degree d the monomial basis functions u“, u%!, ... tend to take 


large values and therefore need to be multiplied with small coefficient values 
Cd, Cd-1,.... This can lead to ill-conditioned matrices, from which numerical 
instability issues can arise [35, pp. 272-274]. Increasing the degree of a poly- 
nomial approximation function also tends to result in an undesired oscillating 
behavior instead of better accuracy [35, pp. 292-294]. This effect is called 
Runge’s phenomenon [150]. 
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z 


Figure 2.2: Curve segment (left) and surface patch (right) [5, pp. 508-509,adapted]. 


The explicit form of a curve y = f(x) in two dimensions, x, y, states the value 
of one dependent variable, y, in terms of the independent variable, x. Most 
curves can also be represented in the implicit form f(x, y) = 0. 


However, curves in three dimensions are not as easily defined in implicit form. 
Instead, they are usually represented in parametric form in terms of an inde- 
pendent variable u, called parameter. An advantage of the parametric form is 
that each dimension is defined by an explicit function that is not coupled to the 
others: 

pu) = [x(u), y), z(u)]" (2.9) 


For example, a cubic polynomial parametric curve p (u) reads 


3 
p(u) = fer =u'c, U= 1, u, u, |; 
k=0 (2.10) 


T T 
Ck = le. Cyk> ca] » C= leo, Cj, C2, cs] 


c is viewed as a 4 x 1 matrix, whose elements cz are 3x 1 matrices. Each 
coefficient vector cg has independent components for each dimension, hence 
there are three independent functions, each of the form (2.8) with d = 3. With- 
out loss of generality, 0 < u < 1 can be assumed. Figure 2.2 illustrates such a 
curve segment. 
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Defining each explicit function in terms of two parameters u and v, leads to a 
parametric surface p(u, v), e.g. for three dimensions: 


p(u, v) = [x(u, v), y(u,v), z(u, v)]" (2.11) 


A surface patch, as also depicted by Figure 2.2, can be defined by varying u 
and v over the rectangle 0 < u,v < 1. A bicubic surface patch is given by 


i=0 j=0 (2.12) 


-fi 2,311 u 
v»=|1,v,v,v > Cij = |Cxij> Cyijo Czij 


G= [ex] is a 4x4 matrix with elements c;;. A surface patch can be seen as the 
limit of a collection of curves that result from keeping one parameter constant 
and varying the other [5, pp. 503-515]. 


Spline functions are piecewise defined, which decouples the number of param- 
eters from the degree of the function. Spline functions allow to increase the 
number of function parameters while still using polynomials of low degree in 
order to avoid numerical stability issues and Runge’s phenomenon. At the join 
points, neighboring pieces of a spline function are continuously differentiable 
to a certain extent. A function that is continuously differentiable up to its r' 
derivative is called a C” continuous function [35, p. 26]. 


The remainder of this subsection reviews the most common cubic spline types. 
Figure 2.3 illustrates them and Table 2.6 compares their features. 


A straight-forward approach to construct a cubic spline function is to connect 
cubic polynomials by specifying continuity requirements at their join points as 
explained in [35, pp. 294-299]. In literature this approach is also referred to as 
natural polynomial spline [9]. 


For example, connecting two cubic polynomials p(u) and q(u) while demand- 
ing continuity up to the second derivative leads to eight degrees of freedom, 
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y  p(t)=a(0) 
p'(1)=q'(0) 
p"(1)=q"(0) 


(ec) Hermite (d) Bézier 
y 
p(1)=q(0) 
p'(1)=q(0) 
p=q  P"(1)=q"(0) 
ETETETT 


(e) B-spline (£) Catmull-Rom spline 


Figure 2.3: Cubic curve segments p (u), and q (u) of various representation types depending on 
parameter u € [0, 1] with their control points p;, qi, i = 0,..., 3 indicated by dots. 
Dotted lines indicate borders of convex hulls defined by the control points. [5, pp. 510- 
534, adapted]. 


which need to fulfill three equality constraints at the join point. Figure 2.3a 
illustrates this case. In order to determine possible spline functions, an un- 
determined linear system of equations needs to be solved. Each additional 
polynomial segment leads to four additional degrees of freedom and three ad- 
ditional continuity constraints. Furthermore, the parameters of all segments 
of the same dimension influence each other via the continuity contraints and 
therefore changing the spline function at any point can affect the shape of the 
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Table 2.6: Cubic representation types for functions, curves and surfaces. 


Continuity at Control point Convex hull 


Representation type Control |. í . f 
join points interpolation property 

Natural polynomial Global Cc? = = 
Interpolating polynomial Local co Yes No 
Hermite Local c! Yes No 
Bézier Local (eu Yes Yes 
B-spline Local Cc? No Yes 
Catmull-Rom spline Local c! Yes No 


function at any other point. This feature is called global control and mostly 
undesired because it requires that all segments are designed within a single 
global calculation. 


In contrast, the spline types described below offer local control of shape, mean- 
ing that changing a parameter only influences the function locally. Therefore 
each segment can be designed individually. Furthermore, they achieve continu- 
ity without solving a system of equations during application. 


Geometrically interpretable parameters called control points are used to design 
the shape of these functions. Assume four control points 


Ppr = [xr Yk, zk], k = 0,1,2,3, (2.13) 


which are equally spaced at u = 0, !/3,?/3,1. c is determined such that p (u) 
from (2.10) interpolates the control points px. The conditions 


Po = p (0), pi = pC), p2 = PCP), p3 = p(l) (2.14) 
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read in matrix form 


po u(0)” 1 0 0 0 
T 2 3 
Rake pe Pı re aa = 1 13 (AB) (/3) 01) 
P2 u (2/3) 12 Ep? P 
P3 ul)" 1 1 1 1 


p and c are 4 x | matrices, whose elements are 3 x 1 matrices. The control 
points given by c = A7'p = Mp with geometry matrix M = A”! achieve 
continuity at join points, but not for the derivatives. Figure 2.3b shows two 
curve segments p(u), q(u) with their control points [5, pp. 509-517]. 


For C? continuity, quintic polynomials would have to be used. Increasing the 
degree of a polynomial segment from three on complicates the calculation 
process with each increment significantly [11]. In particular, for polynomials 
of degree higher than four there is no general closed-form solution [23]. This 
makes nonnegativity checks by determining the roots costly. 


By substituting M into (2.10), p(u) can be expressed in terms of the blending 
polynomials b(u) with 


b(u) =u'M = [bo(u), bı(u), b2(w), b3(w)] (2.16) 
and the control point vector p: 

p(u) =u'c = u' (Mp) = (u'M)p = b(u)p (2.17) 
Interpolating curves can be extended to interpolating surfaces. Assume a 4 x 4 
control point matrix P = [ Pi j| comprising 16 independent three-dimensional 


control points p;;, i =0,...,3, j =0,...,3, which are again equally spaced at 
u, v = 0, 1/3, 2/3, 1. The curve for v = 0 must interpolate poo, P10, P20, P30: 


o 
p(u,0) = u'M | poo, Pio Pa P3) =u" CIL, 0,0, 0)" (2.18) 
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v = 1/3,2/3,1 give the other three interpolating curves. Putting all of them 
together leads to u'MP = u'CAT with A“! = M. Solving for coefficient 
matrix C = MPM" and substituting with (2.16) results in 


3 


3 
p(u,v) =u'MPM'v = b(u)Pb(v)' = > y bi(wb;(w)pij.  Q.19) 
i=0 j=0 


Each term b;(u)b;(v) describes a blending patch. The surface is formed from 16 
blending patches, each weighted by a control point. Surfaces that are created in 
such a way are called tensor product surfaces and are an example of separable 
surfaces. They allow to work with functions in u and v independently [5, 
pp. 510-517]. 


The spline types described in the following can also be represented using (2.17) 
and (2.19). They differ only in the blending functions b(u) from (2.16), because 
their designs fulfill conditions other than (2.14) and (2.18). 


Hermite curves depicted in Figure 2.3c offer C! continuity between different 
segments. However, instead of p1, p2, the user needs to specify the values at 
the join points up to their first derivatives p*(0) and p’(1), so that consecutive 
tangents are collinear. This is often undesired because without an analytic 
formulation of the data, no derivative information is available. 


Bézier curves are approximations to Hermite curves and do not need derivative 
information. They result when the derivatives of two segments at their join 
point are not required to be exactly equal. Therefore cubic Bézier curves only 
have C? continuity. Each segment is within the convex hull given by its control 
points, see Figure 2.3d. For spline types with this convex hull feature, only the 
discrete control points but not the continuous curve itself need to be evaluated 
to ensure that the curve fulfills certain constraints, e.g. nonnegativity. 


For higher continuity, an alternative to increasing the polynomial degree is to 
not require the curve to interpolate any control point. This leads to B-spline 
curves (Figure 2.3e). Each of their segments only spans the distance between 
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its two middle control points. In general they only come close to the control 
points but compensate with C? continuity at the join points. 


Omitting the requirement that each spline segment must lie within the convex 
hull of its control point allows to form other types of splines. One of the most 
popular is the C! continuous Catmull-Rom spline (Figure 2.3f) that interpo- 
lates its control points. Compared to a Hermite curve it does also not require 
derivative information [5, pp. 509-535]. 


The remainder of this work will use cubic B-spline functions. Compared to 
natural splines, they offer numerically stable computations, local control and do 
not require solving a system of equations during their application. Compared 
to the other cubic spline functions, only B-splines are both C? continuous and 
fulfill the convex hull property. C? continuity is an important requirement for 
technical applications such as the definition of a jerk-free trajectory [55]. 


The convex hull property in combination with the geometrically interpretable 
control points offers a convenient and computationally efficient way to enforce 
constraints on the function because no evaluation of the function is needed. 
Enforcing constraints on an approximation function can be beneficial in case of 
knowledge about the data, e.g. that it is nonnegative. 


O Due to their favorable features, B-spline functions, curves and surfaces are 


widely used for data approximation [4, 84, 202] and for defining paths and 
trajectories of vehicles [17, 42, 45], robots [26, 45, 108, 162] and industrial 
machines [69, 200]. Moreover, the B-spline representation is common in com- 
puter graphics [93, 193] as well as in signal processing for filter design [128, 
129, 148] and signal representation [120, 142, 143, 176, 177]. Adopted from 
[78, 80]. 


Methods for B-spline data approximation 


O According to (2.17), the value of a B-spline function equals the sum of basis 
functions (B-splines) or blending functions, weighted with their corresponding 
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control points. Each B-spline in (2.16) is only nonzero within a bounded inter- 
val so that effects of control point changes are only local. Hence, the definition 
range of a B-spline function is also bounded. More extensively these features 
will be described in Section 3.1 and by Figure 3.2. Adopted from [80]. 


o A common application of B-spline functions, curves and surfaces is fitting 
of data points. Fitting can either be interpolation or approximation. An inter- 
polating B-spline function f must pass through all of the data points (sp, Yp), 
Le. f (Sp) = Yp Vp. In contrast, an approximating B-spline function only mini- 
mizes the residuals f (sp) — yp between the function and the data but does not 
pass through the data points in general. 


In offline applications, in which all data points are available at once, fitting 
B-spline functions are often determined by least squares (LS) methods [31, 
113, 159]. With the standard formula in batch form, all data points have to be 
collected first and then processed simultaneously. Therefore the number of data 
points P must be bounded. The computation usually involves a Cholesky or 
QR factorization and requires O(P) operations if one takes advantage of the 
banded matrix structure [20, pp. 327-331]. Such algorithms are stated in [20, 
pp. 117-121] and [58, pp. 152-160]. With the least squares (LS) algorithm each 
data point influences the result to the same extent. The weighted least squares 
(WLS) estimator, stated in (3.14), allows to weight measurements relative to 
each other [113, pp. 119-123]. 


In online applications data points are observed one after another and an ever- 
growing amount of data is common. Two groups of LS algorithms for online 
applications can be distinguished: First, growing memory LS algorithms apply 
an exponential weighting that forgets old data. Second, sliding window LS 
algorithms discard old data completely and need only finite storage [201]. Slid- 
ing window LS and sliding window WLS algorithms are proposed in [83, 201] 
and [33, 184], respectively. Re-computing the fitting function from scratch with 
each new data point is costly. Rank update and rank downdate methods allow to 
reuse an already known factorization for an efficient update of a solution after 
observations have been added or deleted [66, 124, 125]. 
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In online applications data points can be outside the definition range of the 
B-spline function if the magnitude of the data points is not exactly known or 
changes over time, e.g. because one dimension of the data point refers to its 
observation time, that keeps increasing. Data points outside the definition range 
cannot be taken into account for approximation. Thereby the problem arises 
that the approximation might not reflect the data anymore. When using WLS 
the bounded definition range of B-spline functions does not present a problem 
because the number and position of B-splines can be changed if the fitting 
function is re-computed from scratch. Moreover, rank modification methods 
also support adding or deleting matrix columns [66]. This allows to extend, 
shrink or shift the definition range of the B-spline function. Adopted from 
[80]. 


O In nonlinear weighted least squares (NWLS) problems, the solution depends 
on the function parameters in a nonlinear fashion. Based on the results of nu- 
merical experiments, [73] reports that a B-spline function is useful for solving 
NWLS problems as well because of its piecewise polynomial character and 
smoothness. For NWLS problems, several batch methods that work iteratively 
can be found in literature, e.g. the Newton method, Gauss-Newton method, 
Levenberg-Marquardt method, dog leg method of Powell, hybrid method of 
Madsen, Levenberg-Marquardt-Fletcher method. None of the algorithms is an 
exact method that computes an optimal solution [73]. A method for separa- 
ble NWLS problems in which some parameters affect the solution linearly is 
derived in [149]. 


The Levenberg-Marquardt (LM) algorithm, described in more detail in Subsec- 
tion 3.4.1, solves in each iteration a linearized NWLS problem [35, pp. 222- 
224]. A sliding window implementation of the LM algorithm is stated in [38]. 
Adopted from [78]. 


O Recursive methods compute an approximating B-spline function recursively 


meaning that the approximation is updated with each new data point. This ap- 
proach is preferred in online applications, in which data points are observed one 
after another. Recursive algorithms such as recursive least squares (RLS) [164, 
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pp. 84-88] usually require less computational power than batch algorithms be- 
cause they use smaller matrices and vectors whose sizes do not depend on the 
number of data points. The recursive computation is also referred to as progres- 
sive, iterative or sequential [110, 192]. This work will use all mentioned terms 
synonymously. In computer science, however, an iteration differs from a recur- 
sion. In [106], fitting B-spline curves and surfaces are iteratively constructed 
based on the idea of profit and loss modification without solving a linear sys- 
tem. The authors of [36] build on the progressive and iterative approximation 
technique for B-spline curve and surface fitting and prove that their algorithm 
achieves a least squares fit to the data points. A recursive algorithm for optimal 
smoothing B-spline surfaces inspired by the RLS method is presented in [50]. 
Algorithms that involve a Kalman filter (KF) are stated in [67, 104]. Adopted 
from [80]. 


Recursive algorithms for NWLS problems can be based on nonlinear filters. 
Overviews of commonly used filters and their features are provided by [57] and 
[187, p. 615]. In [2] a solution via a modified extended Kalman filter (EKF) is 
investigated and algorithms for offline and online applications are stated. The 
two-step estimator proposed in [68] splits the problem into a linear subproblem 
which is solved with a KF and a nonlinear subproblem to which the Gauss- 
Newton algorithm is applied. Section 2.4 will give an introduction into the 
aforementioned adaptive filters and categorize them. 


Research gap: DI Publications of other researchers regarding the recursive 


data approximation with a B-spline function assume a constant definition range. 
For example, the approaches based on the KF in [67, 104] require that the 
KF state vector contains all control points that are estimated during the whole 
approximation procedure. Therefore the number of control points has to be 
bounded and specified in advance. As a result, these algorithms can only ap- 
proximate data points that are within the bounded B-spline function definition 
range determined at the beginning. Other researchers have not addressed the 
possible issue of data points outside the inital B-spline function definition range. 
Adopted from [80]. 
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Therefore the research gap regarding the recursive data approximation with 
a B-spline function consists of the lack of available online methods for data 
approximation with B-spline functions in the linear WLS and NWLS sense 
that can handle data that leaves the initially chosen bounded B-spline function 
definition range. For example, leaving the initial definition range is needed 
if an additional data point of an unbounded data set is observed in each time 
step and the values in at least one dimension of the data points keep increasing 
because they refer to the observation time. More detailed application examples 
are stated on page 78 and illustrated by Figure 3.9. 


2.4 Adaptive filters 


O As a foundation for Chapter 3 and Chapter 5 this section gives an overview 
of several kinds of adaptive filters with focus on the Bayesian approach to state 
estimation, which calculates the probability density function (PDF) of the un- 
known state of a dynamic system. The required information for calculating the 
PDF stems partly from a system model and partly from previous measurements. 
The state estimation is performed by a recursive filter that alternates between 
a time update that predicts the state via the system model and a measurement 
update that corrects the estimate with the current measurement. Adopted from 
[78]. 


Assume a system whose input u and output y are measurable and known for a 
sequence of time steps p = 1,..., P. 


Adaptive filters include a system model f with free model parameter x. Sequen- 
tially they adjust x based on the sequence of u and y such that the parameter es- 
timate £ minimizes the residuals between the system model output $ = f(&, u) 
and y [111, p. 1]. 


Filtering means that £,, referring to time step p, is calculated based on data of 
the current and all previous time steps, hence p = P. This is also referred to 
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as measurement update. The estimation problem is called prediction in case of 
p > P and smoothing if p < P [72, p. 7]. 


Adaptive filters can be subdivided according to their system assumptions. In 
systems, in which x is assumed to change over time, it is also called system 
state, whereas x is mostly referred to as parameter in systems, in which it is 
assumed to be constant. 


Linear adaptive filters 


o RLS and KF compute an optimal state estimate for systems with linear system 
and measurement equations as well as Gaussian system and measurement noise. 
They differ in that RLS assumes a constant state, whereas the KF is designed 
for tracking a time-variant state [164, p. 129], [57, pp. 3-5]. Use cases include 
parameter estimation [83] and path planning [189], respectively. The KF will 
be explained in Subsection 3.3.2 and applied within RBA in Subsection 3.3.3. 


Nonlinear adaptive filters 


In many scenarios the linear Gaussian assumptions do not apply and subopti- 
mal approximate nonlinear Bayesian filters such as the extended Kalman filter 
(EKF), unscented Kalman filter (UKF) or particle filter (PF) are required [8]. 


The EKF applies a local first order Taylor approximation to the nonlinear system 
and measurement functions via Jacobians, in order to keep the linear state and 
measurement equations. System and measurement noise are both approximated 
with zero-mean Gaussian PDFs [57, p. 52]. Although the EKF is not suitable 
for systems with strong nonlinearity or non-Gaussian noise, it is still often 
successfully used for nonlinear state estimation [27]. For example, NWLS 
approximation via a modified EKF is presented in [2]. 


An alternative to the approximation of nonlinear state and measurement func- 
tions is the approximation of the PDFs. This can be done by propagating few 
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state samples called sigma points through the nonlinear functions. A filter that 
follows this approach is referred to as sigma point Kalman filter. A well known 
representative is the UKF. It uses 2- J +1 deterministically chosen sigma points, 
whereby J denotes the system state dimensions. The PDFs are approximated 
as Gaussians whose means and variances are determined from the propagated 
sigma points [57, pp. 3-5, 62-70]. 


Compared to the EKF, the UKF offers at least second order accuracy [37] and 
is a derivative free filter [57, pp. 62-63], meaning that it does not require the 
evaluation of Jacobians, which is often computationally expensive in the EKF 
[27]. Several publications report nonlinear problems in which the UKF per- 
forms better than the EKF, e.g. for trajectory estimation [37, 63]. However, if 
the PDF cannot be well approximated by a Gaussian, because the PDF is multi- 
modal or has a strong skew, the UKF will also not perform well. Under such 
conditions, sequential Monte Carlo methods like the PF outperform Gaussian 
filters like EKF and UKF [8]. 


The PF approximates the PDF by a large set of randomly chosen state samples 
called particles. The state estimate is a weighted average of the particles. With 
increasing number of particles the PDF approximation by the particles becomes 
equivalent to the functional PDF representation and the estimate converges 
against the optimal estimate [8]. For nonlinear and non-Gaussian systems the 
PF allows to determine various statistical moments, whereas EKF and UKF are 
limited to the approximation of the first two moments [27]. However, the num- 
ber of particles that is needed for sufficient approximation of the PDF increases 
exponentially with the state dimension [122]. The PF has been successfully 
applied to optimization [99] and prediction [190] of trajectories as well. 


Many use cases involve a mixed linear/nonlinear system. Typically there are 
few nonlinear state dimensions and comparatively many linear Gaussian state 
dimensions. The marginalized particle filter (MPF) is beneficial for such prob- 
lems as it combines KF and PF. The PF is only applied to the nonlinear states 
because the linear part of the state vector is marginalized out and optimally 
filtered with the KF. 
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This state decomposition is called Rao-Blackwellization [72, pp. 49-50] and 
can be described as an optimal Gaussian mixture approximation. Therefore the 
MPF is also called Rao-Blackwellized particle filter [70]. Marginalizing out 
linear states from the PF strongly reduces the computational effort because less 
particles suffice. This often enables real-time applications. Simultaneously the 
estimation accuracy usually increases [27, 205]. When the state dimension is 
large, the MPF tends to outperform the PF [123]. 


In the recent past, several publications have proposed approaches for localiza- 
tion [135, 191] and trajectory tracking [109, 205] that are based on the MPF 
because of its advantages for mixed linear/nonlinear systems. Automotive use 
cases include a road target tracking application, whose multi-modality requires 
using a PF or MPF [166]. The MPF is chosen as it allows reducing the number 
of particles for less computational effort. Similarly, [122] presents aMPF appli- 
cation for lane tracking, in which the achieved particle reduction compared to 
a pure PF enables executing the algorithm in real-time in an embedded system. 
Adopted from [78]. 


According to (2.17), between spline function value and spline function control 
points there is a known linear relationship given by the blending or basis func- 
tions. Section 3.1 will state in (3.3) and (3.5) that this applies to both the value 
of a B-spline function and its derivatives as well. Therefore NWLS approxima- 
tion leads to a mixed linear/nonlinear problem as long as there are target criteria 
that refer to the B-spline function or its derivatives directly. For being able to 
take into account this known linear relationship instead of having to estimate 
it, the iterative algorithm for NWLS approximation, NRBA, that is defined in 
Subsection 3.4.3, includes an MPF, which is described in Subsection 3.4.2. 


Kernel adaptive filters 
The two previous subsections assumed that knowledge about the system state 


with respect to the measured input u and measured output y is available. Some- 
times this does not apply or the system features are difficult to specify. In such 
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cases the unknown system function f with y = f(x) can be approximated by f. 
The approximation f is also called black-box model because it is created only 
from the data itself without system knowledge. 


Ôf can be determined with a kernel adaptive filter (KAF). Using a nonlinear 
mapping ®(-), KAFs transform the input data u into a high-dimensional feature 
space, in which more degrees of freedom are available to solve the problem 
with a linear adaptive filter. Due to a property of reproducing kernel Hilbert 
spaces no costly explicit mapping is needed to compute inner products (+, -) 
of the transformed data ®(x;). Instead, inner product algorithms can be per- 
formed implicitly in feature space by replacing all inner products in the original 
problem by a kernel function x [146]: 


K(X j, Xm) = (D(X ), P(Xm)) (2.20) 


Often Gaussian kernels are used. These kernels differ from multi-dimensional 
PDFs of the standard normal distribution with variance o only in the missing 
normalization constant 1/ (o V27): 


(xj, Xm) = exp ||x; - xml] 120%) (2.21) 


The nonlinear mapping f is calculated as a linear combination of kernels, each 
of which is weighted by its coefficient. For example, by transforming the RLS 
cost function into feature space, a nonlinear kernel RLS (KRLS) algorithm is 
derived. With most kernel-based methods the required memory increases with 
the processed data points. In contrast, the Fixed-Budget KRLS (FB-KRLS) 
algorithm can operate with constant memory by pruning the least significant 
data [180]. This enables long-term online applications. 


In summary KAFs are nonlinear adaptive filters that result from applying kernel 
method concepts to linear adaptive filters. KAFs combine the capability of 
approximating any nonlinear function with the advantage that they only require 
solving a convex optimization problem instead of a nonlinear problem. For 


37 


2 Scientific and technical state of the art 


convex problems there are solvers that operate within a fixed time interval. This 


is important for real-time applications [146]. 


Subsection 4.7.2 will apply the FB-KRLS to represent the power train of the 
research vehicle. The model is updated constantly during vehicle operation. 
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3 Data approximation with 
B-spline functions 


This chapter deals with the general problem of recursive least squares approx- 
imation of data points with a B-spline function. First, Section 3.1 states a 
B-spline function definition in matrix form. Section 3.2 introduces a data set 
which will be approximated by a B-spline function. Section 3.3 addresses 
the case of linear least squares approximation and Section 3.4 the nonlinear 
case. Each of the two latter sections proposes a novel recursive algorithm 
and compares it with a well-known method for batch processing. Section 3.5 
summarizes the scientific contribution. 


3.1 B-spline function definition 


A B-spline function consists of several polynomial basis splines (B-splines), all 
of which have the same degree. A B-spline of degree d is piecewise defined 
using d + 1 polynomial functions of degree d. Figure 3.1 illustrates the compo- 
sition of a cubic B-spline. The solidly drawn parts of the depicted polynomial 
functions pı, p2, p3 and pa build a B-spline. 


Let [k, K2] denote a closed interval between and including x; and x2 and let 
(K1, K2) denote an open interval between but not including x; and x2. 


In the interval [x1, «2), the B-spline is given by the first polynomial pı, in the 
interval [k2, K3) it is given by the second polynomial pz and so on. Outside the 
interval (kı, k5) the B-spline is zero. The dashed curves are only shown for 
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Figure 3.1: Piecewise definition of a cubic B-spline by the solid parts of four cubic polynomials 
Pi» P2, - - -, p4. Vertical lines indicate the knots kj, ..., K5. 


better visualization of the polynomials but are not part of the B-spline. The 
interval borders K1, K2,..., K5 are also called knots. Their distances influence 
the shape of the polynomials and therefore the shape of the B-spline as well. 


O A B-spline function is also piecewise defined. Its value is given by the 
weighted sum of J B-splines of degree d. k with K = (K1,K2,..., KJ+d+1) iS 
the knot vector. Strictly increasing knot values (kx < Kk+1, k = 1,2,..., J +d) 
are assumed. x and d determine the number and shape of B-splines. The j-th 
B-spline b;(s), j = 1,2,...,J is positive only for s € (Kj, Kj+a+1) and zero 
elsewhere [113, pp. 37-42]. 


The following definitions originate from [113, pp. 47-50, 65-70]: Let [Ky Ku+1) 
be a spline interval and let u denote the spline interval index with d+1 < u < J. 
For s € [Ku Ku+1), the B-splines b;(s), j = u -d,..., u can be nonzero. Their 
values for a specific s € [Ku, Ku+1) can be summarized in the B-spline vector 
bu,a(s) = (bu-a(s), bu-a+1 (8), . . -,bu(s)) € R!*(4+) which can be computed 
according to (3.1): 


bua(s) = By i(s) Bu2(s)... Bus (8)... Bua(s) 
ER!x2 ER?X3 EROX(6+1) eRAX(d+1) 
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The B-spline matrix By,5(s) € R°**) is defined for each ô € N with ô < d 
and given by 


B5(5) = 
Ku+1=8 S—Ky+1-6 
Ku+lKu+l-5 Ky+1—Ky+1-5 0 == g 
Ky+2—S SK 42-6 
0 Kyut2—Kyt+2-5  Ku+27Ku+2-5 277 0 (3.2) 
0 0 Ku+ó—58 S=Ku 


Ku+ó "Ku Ku+s "Ku 


The B-spline function f : D > R,s» f(s) has the definition range D = 
[K4+1, KJ+1). For s € [Ku Ku+1), the B-spline function is given by 


FCs) = dy. a(S)X p,a (3.3) 
with control point vector 
= 
Spa = (xid ER E (3.4) 


Adopted from [78, 80]. As stated, the spline types mentioned in Section 2.3 
can all be defined using the same structure because they only differ in the 
geometry matrix M. Therefore (3.1) and (3.3) correspond to (2.16) and (2.17), 
which were obtained for an interpolating polynomial curve. 


Figure 3.2 illustrates the construction of a B-spline function. Each cubic 
(d=3) B-spline b; is weighted with a corresponding control point x;, here 
xı = 1.5, x4 = 0.5 and all other control points equal one. The cubic B-spline 
function f(s) is given by the weighted sum of B-splines but only defined in the 
interval [x4, K7). Its shape is given by the solid part of the black curve. 


The black dots denote the control points of the B-spline function. The horizontal 
position of the j-th point is determined by the value of s, for which the j-th 
B-spline reaches its maximum. The vertical position is identical to x;. In each 
interval [kx, Kx+1), k = d+1,...,K, the B-spline function lies within the convex 
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~~) A 


definition range of f(s) 


Figure 3.2: Construction of a cubic B-spline function: Equidistant knots «1, K2, ..., K10 (indicated 
by vertical straight lines) and J = 6 resulting B-splines b;, j = 1,2,..., J weighted 
with control points x1, x2, ..., X6. The black line indicates the sum of all weighted 


B-splines and its solid part is the B-spline function. Black dots represent the control 
points and the gray area is the convex hull formed by the first four control points. 
Adopted from [80]. 


hull of the relevant control points. For s € [k4, Ks), the four leftmost control 
points are the relevant ones. They form the convex hull indicated by the gray 
shaded area. 


o A B-spline function of degree d is d — 1 times continuously differentiable. 
o” 
Os” 


For r € No, the r-th derivative 25 f(s) of the B-spline function with respect to 


s is given by 


o” or 
sr f(s) = gyr nal IX pd (3.5) 
with B-spline vector 
o” 
Bar teal) = 
ae Bus). Baar Bl ar Bip ifr <d (3.6) 
Oix(d+1) otherwise. 
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3.2 Structure of the data set 


O1x(a+1) denotes a 1 x (d + 1) zero matrix. The matrix Bi as ROXO+) jg 
obtained by differentiating all entries in B,,s(s) with respect to s: 


-1 1 0 
Ku+lKu+l-58  Ku+lTKu+l-ó Be 
, . F 4 é 
Bis = E e ie ; (3.7) 
0 -1 1 


Ku+ó “Ku Ku+ó “Ku 
Adopted from [78, 80]. 


According to [161] the antiderivative or indefinite integral F(s) of f(s) with 
degree d'™ = d + 1 is given by 


F(s) = f f(s)ds = bim gm kS Xq (3.8) 
b™ uses the stricly increasing knot vector x! which differs from « only by an 
additional knot at either end for the increased degree. u™ equals u + 1. The 


control point vector x! comprises one additional component than x and can 
be computed according to (3.9): 


J 
K54 gin — Kj 
Int _ Int _ tal J 3 an 
Xi =0, de ES TOTES EN pe O | (3.9) 
i=1 


The definite integral of f(s) over the interval [s¡, s2] then reads 


1 f(s)ds = F(s2) — F (s1). (3.10) 


3.2 Structure of the data set 


O Section 3.3 and Section 3.4 will work with the set {(Sp.¥p)}p=1,2,....P of P 
data points. p denotes the time step, at which data point (sp, Y p) is measured 
or observed. 
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Sp is the value of the independent variable s at time step p. The vector of 
independent variables s is given by s = ($1,...,5p. -> sp)”. 


Jp = IS A JAN is a vector of V, measurements y that 
refer to s, and may come from different sensors. V, € N can be different for 
each y, but it is assumed that V, << P Yp. The vector of all measurements y is 
composed as follows: 


Ps VE VS) (3.11) 
tee E ee ————- 


=y] =y p 


Adopted from [78, 80]. 


3.3 Methods for linear weighted least 
squares problems 


Subsection 3.3.1 describes the WLS approach followed by the KF algorithm in 
Subsection 3.3.2. Subsection 3.3.3 presents the RBA algorithm. Its effective- 
ness is demonstrated in comparison with the WLS solution in Subsection 3.3.4. 


3.3.1 Weighted least squares estimator 


O The linear weighted least squares (WLS) method estimates the constant state 
vector x € REX! of a linear system 


y =Cx+v. (3.12) 


ye RAX! is the vector of measurements from (3.11) and C denotes the mea- 


surement matrix that relates x to y. The measurement noise v € RN*! is 
assumed to be an uncorrelated white noise process with mean zero. This im- 
plies that the covariance matrix of measurement noise R is a diagonal matrix 


and R;.;, i = 1,..., N is the variance of measurement y, which can differ from 
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the variances of other measurements [164]. The assumptions for {v} can be 
generalized to a correlated noise process. This is termed generalized linear 
model [138, p. 143]. Then R is a positive definite matrix [20, p. 374]. 


The linear WLS estimate X minimizes the sum of squared errors between the 
measurements y and the vector Cx which are weighted with the reciprocals of 
the variances of the measurements: 


% = arg min(y - Cx) R! (y - Cx) (3.13) 


The solution to optimization problem (3.13) is given by the closed-form estima- 
tor 
ê = (CTR IC) ICTR !y [164]. (3.14) 


From (3.3) follows that the value of a B-spline function is a linear combination 
of its control points. Therefore WLS can be used to determine the control 
points such that the function approximates the set of data points defined in 
Section 3.2. Then C is a (S Vp) x J matrix because y comprises pe Vp 
scalar components yp y, P = 1,..., P, v=1,..., Vp (c.f. (3.11)) and there are 
J B-splines. yp,y is the v-th component of y (5 = a Va + v) and provides 


information about = f (Sp) with sp € [Ku Ku+1) and anr € No. The P-th row 
of C is given by C;:1....y = c with 
or 
c = |Oixu-(d+1)}» gyr Pea (Sp)s Dix) (3.15) 


Herein 0, denotes ar x c zero matrix. Adopted from [80]. 


3.3.2 Kalman filter 


O The Kalman filter (KF) is an established method for estimating the state of a 
dynamic system. Applications include tracking, navigation, sensor data fusion 
and process control [60, pp. 4-5]. The KF can be seen as a generalization of the 
RLS method [164, p. 129]. 
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The linear KF estimates the state vector x, € RX of a linear time-discrete 


system 
Xp = ApXp-1 +B, 4, + wp (State equation) (3.16) 
Yp =CpXp+Up (Measurement equation) (3.17) 


where p € N denotes the time step. Ap is the state transition matrix that relates 
Xp-1 t0 Xp, Up € R™*' is an input signal vector with known influence on x p 
and B, is the input matrix that relates up to xp. The vector of measurements 
is denoted by y,, € RN*! and C p is the measurement matrix that relates x, to 
Yp- Op € R£*! is the process noise with covariance matrix Q, and v, € RN*! 
is the measurement noise with covariance matrix Rp. Both {wp} and {vp} are 
uncorrelated white noise processes with mean zero which implies that Q, and 
Rp are diagonal matrices [164, p. 124]. 


The KF consists of a sequel of equations, which are computed for each time 
step and summarized in Algorithm 1, in which J denotes the identity matrix 
with appropriate dimensions. The KF performs a time update followed by a 
measurement update. During the time update, the state estimate is updated 
based on the knowledge about the system specified by (3.16). Both the a priori 
estimate X,, and the covariance P, of the a priori estimation error are calculated. 
During the measurement update, the Kalman gain X, is computed and used 
together with the information provided by measurement y,, for the calculation 
of the corresponding a posteriori quantities iF and P? [164, pp. 124-129]. 
KF generalizations for correlated or colored noise processes are stated in [164, 
pp. 183-193]. 


If the state vector x, is constant, then Ap = I, wp = 0 and up = 0, whereby 0 
is the zero matrix with appropriate dimensions. In this case, the time update is 
redundant and the KF simplifies to the RLS algorithm [164, p. 129]. 


Data approximation by a polynomial using RLS is described in [164, pp. 92-93]. 
Due to the unbounded definition range of a polynomial, applying RLS saves 
computational effort without limiting the approximation compared to a KF. 
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Algorithm 1: Kalman filter. Adopted from [80]. 
Input: £7 Php Up, Yp Ap, Bp, Cp, Qp Rp 
/* Time update x / 
ee Ap? 1 +B,u, 
2 P, APP ¡Ap' +Qp 
/* Measurement update x / 
3 Kp  P5Cp (CpP5Cp' + Rp)" 
4 £5 = 85+ Kp(Vp -Cp3,) 
Pi E (I -KpCpP E -KpCp) +KpRpKp" 
Output: $$, P} 


= 


u 


Data approximation by a B-spline based on the KF is proposed in [67, 104]. 
These approaches require that the state vector x, comprises all control points 
that are estimated during the whole approximation procedure. This means that 
Xp, i.e. what is estimated, is constant and only the estimation X, can change. 
Hence, RLS would suffice. Combined with the bounded definition range of a B- 
spline function, the missing ability to change which control points are estimated 
leads to a definition range of the approximation that is also bounded, constant 
and needs to be specified in advance. 


In contrast to these approaches, the recursive B-spline approximation algorithm 
proposed in Subsection 3.3.3 takes advantage of the time update, that the KF 
provides, to shift the control points in xp. As a consequence, RBA can shift the 
definition range of the B-spline function to consider data points outside of the 
current definition range. Adopted from [80]. 


3.3.3 Recursive B-spline approximation algorithm 


O The recursive B-spline approximation (RBA) algorithm computes an approx- 
imating B-spline function f(s) of degree d for the set of data points from 
Section 3.2 iteratively using the KF. Algorithm 2 summarizes the calculations. 
I € N denotes the constant number of spline intervals of f(s). The KF state 
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e e e olo/o|o|o ojofojojo 
e ojofojo ojofojojo ofolojojo ojofojojo 
e ojofojo ojofojojo olo/o|o|o ojofojojo 
. ? ojofojo >: ojofojojo > Tololo lolo > Cololofolo 
e ojofojo ojofojojo olo/[o|o|o ojofojojo 
e e ojofojojo e 
(a) ©) (c) (d) (e) 


Figure 3.3: Changes of covariance matrix elements for d = 3 and I = 3. e indicates a large 
positive value, empty cells indicate zeros and o denotes comparatively small values. 


Adopted from [80]. 
estimate ĉp = (£p,,£p,,--.,£p,)" comprises J = d + I components which are 
the estimated control points of f(s). The knot vector Kp = (Kp; Kp» +++» Kpg) 


for time step p has to contain K = J + d + 1 knots. Dp = [Kp ¿,¡> Kpy,,) 18 the 
definition range of f(s) at p. 


Initialization 


X, is initialized with x = x1yz¡, where 1y.¡ denotes a J x 1 matrix of ones 
and x a scalar quantity of the magnitude of measurements y,,, that refer to 
Z f(s) with r = 0. 


The covariance matrix of a posteriori estimation error P* is initialized with 
Po = PI ¡x3, where Ij, 7 is a J x J identity matrix. The scalar p should be 
chosen large (e. g. 10*) because then X, will quickly deviate from its initial 
value £5 in such a way that f(s) approximates the data points. If the elements 
in P, are small, this signals to the KF that the state estimate £, is very certain 
and therefore it will hardly be updated using the measurements. If the KF 
updates 7, as intended, the elements in Po become smaller as p increases. 


In the long-run, P, is strongly influenced by the covariance matrix of process 
noise Q, according to line 2 of Algorithm 1. If the elements in Q, are large, 
the elements in P, remain large too. This can lead to volatile state estimates 
X» that do not converge to a certain value. Then f(s) will not approximate the 
data points well. For that reason a very small positive value q (e.g. q= 10712) 
is chosen for Qp = GI yx. 
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Algorithm 2: Recursive B-spline approximation. Adopted from [80]. 


vsonaurun u 


N DD DD DREI Fe EF EF EF SE j 
Sa wu Y Ne Oo! FA HAN BUND KF S 
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Input: Kp-1> De P_i Rp, Sp» Jp Kp, X, P, q 
Je length(+7_,) 
K — length(x,-1) 


deK-=J-1 
I-J-d 
o-0 


if sp > Kp-1,,, then 
if sp > Kp-1, then 

| oa-d+l 
else 

| o such that sy € [ 
end 
else if sp < kp-1,,, then 
if sp < Kp-1, then 

| ao -(d+ 1) 
else 

| o such that sp € [Kp-1 41,5 KP- ) 
end 


A Kp-lastiase) 


end 

Ap € R” from (3.19) 

Qp aL 1x3 

ifo > 0 then 

Kp from (3.18) 

Up — (Oix(J-0) xo)" 
O CB M=J-0+1,J-0+2,...,J 
else 

Kp from (3.25) 

Up E (lixo) Dix(s+0)) * 
Qpmm Ps m= 1,2,..., -0 


end 

u such that sp € [Kp > Kp,,,1) 

Vp = length(y„) 

Cp € RYP” from (3.27) 

[8 P5] — Algorithm Le Pe ie Up. Y p Ap, 3x5, Cp, Qp, Rp) 
Output: Kp, +5, P3 
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Time update with shift operation 


RBA compares the knot vector K,-¡ with sp in order to determine whether 
Sp € Dp-1. If necessary, a shift of D,_ı is performed during the time update 
of the KF such that sp € Dp. The variable o indicates the shift direction of 
D,-ı and the number of positions by which elements are shifted. 7 > 0 means 
a right shift of D,-1, © < 0 a left shift of D,-ı and for a = O no shift is 
performed because sp € Dp-1. Algorithm 2 computes o from line 5 to line 18. 


For example, assume d = 3, J = 3 and kp-1 = (1,2,...,10), then D,-1 = 
[4,7). If sp = 8.5, two additional knots are needed to be able to perform 
a right shift by two elements (a = 2). 11 and 12 as additional knots give 
Kp = (3,4,...,12) and hence sp € Dp = [6,9). 


Algorithm 2 distinguishes between a > O ando < 0. It assumes that 
the |o| additional knots are the o last components of the knot vector Kp = 
(Kp » Kp ) « - => Kpg ) in case of o > 0 and that they are the —c first components 
of Kp in case of o < 0. 


Case 1: Right shift of definition range (c > 0) 


The new knot vector is 


Kp E (Kpt oo Kp pao ++ +9 Kp 1g 


(3.18) 


Kpx-o41 KPK-o42? °°? Kpg). 


The elements in 1 are shifted by line 1 of Algorithm 1 using the state 
transition matrix 


l, ifh=g+ 
Ap ER”! with Ap, = pay (3.19) 
j O, otherwise. 
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X, = Apts) updates the old estimate ry to 


A A A A T 

Xp = (CS Xp-1 547 E Xp-1,4>01x0) © (3.20) 
With the second part of the instruction (X, — £, + Bpup), input matrix 
Bp = I ¡xy and input signal vector 


Up = (Mixo) Fixe)" (3.21) 


arbitrary initial estimates x in X,, can be obtained: 


A 


me x = T 
lr ri) (3.22) 
Pes is updated during the time update as well by line 2 of Algorithm 1. The 
first part of the instruction (P, — ApP), Ap") leads to Ph, yo yg T 
p- 


P- Li oes O i 


2,...,J of Pp equal zero. Especially zeros on the main diagonal prevent that 


and all elements in the rows or columns J- o +1,J-0+ 


Ep. pr. > Xp, become different from the initial value x. 


These zeros in P, can be replaced with large values using the second part of 
the instruction, P, — P; + Qp, with 


P ifh=gaQ 


0, otherwise 


Qp ER” with Q 


Peh 


with 
h>2J-0+1, ifo>0 
Q= (3.24) 
h<-o, ifø <0. 
Figure 3.3 depicts different states of P* and P”, respectively. d = 3 and I = 3 
lead to a 6x6 matrix P*. The diagonal values of Pù are initialized with a large 
positive value as depicted in (a). All other elements of Pf are initialized with 
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zeros. After some data points in the second spline interval have been processed, 
different comparatively small values are in the submatrix P3.5 „., (b). 


After data points that fall into the third spline interval have been taken into 
account, only the elements in the first row and column of P* still have their 
initialization values (c). If only the first part of the update instruction is executed 
during a right shift by one element (o = 1), the elements in the last row and 
column of P~ become zero (d). With the second part of the instruction, these 
elements can be set to nonzero values. For Q¢,6 = e and all other elements of 
Q equal zero, matrix (e) is obtained. 


Case 2: Left shift of definition range (c < 0) 


The new knot vector is 


Kp = (epi: Kpo ++ > Kp_ > Kp-1;> Kp-1p +0 > ké) (3.25) 


and the input signal vector reads 


up = Chica bern) : (3.26) 


Effect of shift operation 


The shift operation is the distinguishing feature of RBA compared to the al- 
gorithms in [67, 104]. Due to the shift operation, the total number of control 
points that will be estimated during the application of RBA does neither have 
to be known in advance nor to be bounded. Using Kp and Lo the determined 
B-spline function f(s) can be evaluated at any s € [Kpa ¡> Kp4,7,1)- 

+ 


p-\? 
constant. In Algorithm 2 |o] is chosen just large enough that the current data 


By shifting the entries in Kp-1, 5 and P the required memory is held 


point can be taken into account during the measurement update. The reason 
for this is that the shift operation comes at the cost that parts of an already 
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computed approximation are forgotten unless the values of Kp-1, 51 and 
Posi are stored separately from Algorithm 2 before they are removed from the 
vectors or matrix, respectively. 


Measurement update 


During the measurement update, the information provided by (sp, y p) is used 
to update f(s). With the covariance matrix of measurement noise Rp, different 
components of the measurement vector y,, can be weighted relatively to each 
other. Rp € RV? XV? is a diagonal matrix with positive elements on its diagonal. 
The smaller an entry Rp,,., is, the greater is the effect of the v-th component of 
the measurement error (y, — Cp*,) on Lp: 


The measurement matrix C, is a Vp x J matrix. yp, is the v-th component 
of y, and provides information about Zf (5p) with sp € [Ku Ku+1) and an 
r € No. The v-th row of Cp is given by 


Co... = € (3.27) 


with c from (3.15). According to Cp, (Sp, y p) influences only the estimates 
x Pu-a’ Xp, ig Xp, However, other estimates can still be updated by the 
KF using the information stored in P,,. Adopted from [80]. 


3.3.4 Effectiveness of recursive B-spline 
approximation 


O A numerical experiment in [80] demonstrates the effectiveness of Algorithm 2 


in comparison with the corresponding WLS solution. 
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The data set {(Sp, Y p))p=1,2,..., P from Section 3.2 with P = 5000 is chosen, 
whereby 


Sp = 0.01 + 0.02(p — 1), (3.28) 
20, if 30 < sp < 70 
Yp,1 = and (3.29) 
10, otherwise 
Yp2 = Yp,3 = 0 Yp. (3.30) 


A cubic (d = 3) B-spline function f(s) approximates the data set, whereby it is 

assumed that the measurements y, 9 refer to the slope 2 f(s) of the B-spline 
2 

function and the measurements yp,3 to the curvature 5 f(s). The reciprocals 

of the relative weights between the different target criteria are specified by the 

pia = l, Rove = 

107? and Ro; = 107°. The diagonal measurement covariance matrix of WLS 


diagonal measurement covariance matrix R, € RAS with R 


R e R3?*3? has analogous values: 


-2 3 
Rias = 1, Rigo = 1005 Ri4z;43 = 10 >, 


(3.31) 
i=3(r-1),r=1,2,...,P 


The chosen weighting helps to prevent overshoots and oscillations of f(s) and 
leads to a B-spline function that smooths the jumps of y,,ı in the data set. The 
parameter values are g = 10712, & = 0 and p = 10%. 


In order to investigate the effect of the choice of 7, four runs of RBA with 
I = 1,3,7 and 20, respectively, are performed. kg = (-15, 10, 15,20) is used 
for J = 1, ko = (-15,10,...,30) for I = 3, Ko = (-15,10,...,50) for I = 7, 
and Kg = (15, 10,..., 115) for J = 20. For J = 20 the resulting Do comprises 
all s, of the data set and therefore no shift operation is needed. For J = 1,3 
and 7, RBA has to perform several right shifts by one element in order to be 
able to process all data points. For each shift operation, an additional knot kp 
has to be defined in the vector Kp. Kpg = 25,30,..., 115 is chosen for J = 1, 
Kpg = 35,40,...,115 for I = 3 and kp, = 55,60,...,115 for J = 7. For 
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u (Sp, Yp,1) — WLS 


— RBA;-3 --- RBA;-7 
y2.11f() 


Figure 3.4: Effectiveness of RBA: 40 of the 5000 data points (sp, yp,1) (black dots) and knots 
0, 5, ..., 100 (vertical dashed lines) are shown as well as the B-spline function f(s) 
determined by WLS and RBA with J = 1, 3, 7 and 20. Arrows visualize the definition 
range of f(s) while data points in the interval [95, 100) are processed. Adopted from 
[80]. 


evaluation purposes X,-ı, and Kp-1, are stored separately from RBA before a 
shift operation is performed. 


Figure 3.4 displays the results. As both the data set and knot vector are symmet- 
rical to a vertical straight line through s = 50, the WLS solution is symmetrical 
as well. The RBA solution converges to the WLS solution as / increases. For 
I = 1 and / = 3, the resulting B-spline function is asymmetrical with respect 
to a straight vertical line through s = 50. For J = 7, RBA provides almost 
the same result as for Z = 20. Consequently, in this example J can be reduced 
from 20 to 7 without noticeably worsening the quality of the approximation. 
Lowering / leads to less computational effort in the KF because P,, P and 
Q, are (d +1) x (d +1) matrices and therefore the asymptotic time complexity 
of each individual iteration of Algorithm 1 is O((d +I )3) if the standard method 
for matrix multiplication is used [25, 97]. Under the same conditions both RLS 
and the known methods of other researchers based on the KF need 20 x 20 
matrices because shift operations are not possible. 


Due to the large possible savings in computational effort, RBA is also beneficial 
in offline applications with known finite data set and therefore a definition range 
that is known in advance. 
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In further experiments in [80] Ko is chosen such that the first data points lie in the 
rightmost spline interval and runs with J = 3, J = 7 and J = 20 are performed. 
The deviations of the resulting control point vectors from the those derived 
when the first data points are in the leftmost interval are close to the machine 
accuracy indicating that the effect of the shift operation on the approximation 
result is negligible. Adopted from [80]. 


3.4 Methods for nonlinear weighted least 
squares problems 


Subsection 3.4.1 describes the Levenberg-Marquardt (LM) algorithm followed 
by the marginalized particle filter in Subsection 3.4.2. Subsection 3.4.3 presents 
the nonlinear recursive B-spline approximation algorithm for nonlinear weight- 
ed least squares problems. Its effectiveness is demonstrated in comparison with 
the LM solution in Subsection 3.4.4. 


3.4.1 Levenberg-Marquardt algorithm 


Consider the problem 
& =argmine'e (3.32) 
x 


with the error function e given by 
e = VR”! - (y - ġ(s, x)). (3.33) 


s denotes the vector of independent variables, y the vector of measurements 
and R the covariance matrix of measurement noise. In (3.33) the difference 
between y and & is weighted with the square root of the reciprocal of the 
covariance matrix of measurement noise R. @ is a function whose values 
depend nonlinearly on its control points which are summarized in the control 
point vector x. Due to this nonlinear relationship (3.32) is a NWLS problem. 
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The Levenberg-Marquardt (LM) algorithm solves this problem iteratively. %! 
denotes the known approximation of the solution £ of (3.32) in iteration i. LM 
approximates the error function linearly using 


e(&) zer) +e (I — 8) (3.34) 
and updates £’ according to 
(3.35) 


with correction step size 6. The optimal step size 6 is derived by solving the 
linear LS problem 


3" = arg minlle’(#')6 + eD? + Allo" 2 (3.36) 
Si 


which includes a damping parameter A > 0. The solution to (3.36) with 
ni A fi =f E ® 
Ò = - [ea Tea) + V7] RNA) (3.37) 


is determined by setting the first derivative of its optimization function with 
respect to 6 to zero. A too large 6 can be avoided by choosing A appropriately 


because of ; 
i le^ |l 
lô" ll < o (3.38) 
The LM algorithm terminates if 
[ler (20 Teh < e (3.39) 


with a specified tolerance e. For convergence A needs to be sufficiently large. 
However, a large A results in a small correction step size and slow convergence 
if the solution is still far away from the optimum. Therefore A is controlled 
using heuristic criteria [35, pp. 222-224]. 
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The computational effort of a single LM iteration roughly equals the effort of 
the WLS method. The total computational effort of LM depends on the number 
of performed iterations and is given by O(e?) [175]. 


3.4.2 Marginalized particle filter 


O The marginalized particle filter (MPF) is an iterative algorithm for estimating 
the unknown state vector x, of a system at time step p € N. 

; Peia T my” 
In the MPF, xp is subdivided into x, = (=) , (xX) ) , whereby the KF 
optimally estimates the linear state vector x and a PF estimates the nonlinear 
state vector xy . Exploiting linear substructures allows for better estimates and 
a reduction of the computational effort. Therefore, the MPF is beneficial for 
mixed linear/nonlinear state-space models [155]. Due to Equations (3.3) and 
(3.5), linear substructures will occur in approximation problems as long as there 
are target criteria that refer to the value of the B-spline function or its derivatives 
directly. 


MPF algorithms for several state-space models can be found in [155] along 
with a MATLAB example that can be downloaded from [156]. An equivalent 
but new formulation of the MPF that allows for reused, efficient, and well- 
studied implementations of standard filtering components is stated in [70]. 


For a NWLS approximation, the following state-space model derived from [70] 


is applied: 
EN 41 = AY xy + wy + un (Nonlinear state equation) (3.40) 
4 = Are + ws + už (Linear state equation) (3.41) 
yp = C x +c (xy ) +Up (Measurement equation) (3.42) 


L and Y indicate that the corresponding quantity refers to 


The superscripts 
linear or nonlinear state variables, respectively. A, denotes the state transition 


matrix, u, is the known input vector, y is the vector of measurements, Cp 
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is the measurement matrix, and c is the nonlinear measurement function that 


N 
P 


covariance matrix Q5, wn is the process noise of the nonlinear state vector with 


depends on x}. ws, denotes the process noise of the linear state vector with a 


a covariance matrix QN, and vp is the measurement noise with a covariance 
matrix Rp. The model of the conditionally linear subsystem in the KF has the 


yy T 
state vector (e, (+2) ) , whereby £ describes the linear dynamics of x: 


N N N 

pri = 0 A, Ep i Up + Vp 
L L L L L 
X51 0 A u, wp 


(3.43) 


yp = (o cp) Sp +0 (xp) +0) 
Xp 


The covariance matrix of process noise is and 0 denotes a zero 


L > 
p 
matrix with a suitable size. 


A PF with the model 
N - N 
x =O 
Pere. oe (3.44) 
Vp = Up 


deals with the remaining nonlinear effects. The noise depends on the estimates 
indicated by ^ from the conditionally linear model: 


3.45 
p ~N (c (x) +C, (1) 24", 5,) vn 

with 
Sp =C ¿PF C +Rp (3.46) 


where the superscript 
time update, which is based on the state of (3.40) and (3.41). In contrast, * 


refers to a priori quantities that are computed in the 
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denotes a posteriori quantities that are calculated in the following measurement 
update based on the measurement of (3.42). 


Per and PA are the covariance matrices of the estimation errors that belong 
AL E s 
to X, and & „, respectively. 


The PF uses multiple state estimates, called particles, simultaneously. The 
superscript ? with p = 1,..., P is the particle index and P is the particle count. 
In general, a KF is used for each particle. In the chosen state-space model, 
however, AF ; AY 5 Qs, and o are independent of Yo and x . This implies 
that Re and PS are identical for all KFs, which reduces the computational 
effort substantially [70, 155]. 


Algorithm 3 states the equations for one MPF iteration and was derived from 
[70, 155]. For an implementation in MATLAB, the example from [156] was 
adapted. Note that in Algorithm 3 the measurement update of the previous time 
step p — 1 occurs before the time update for the current time step p, similar to 
the algorithm in [6]. 


In line 4 of Algorithm 3, linear particles are resampled according to their cor- 
responding normalized importance weights. After resampling, particles with 
a low measurement error occur more frequently in the set of particles. Subse- 
. AL, +,p 
quently, all particles 8,1 


by calculating their mean. 


are aggregated in line 5 to a single estimate 1 


After both KF and PF have been time updated, the KF is adjusted based on 
the PF estimates in a mixing step with the cross-covariances of the estimation 
errors, Pe and PS. 


In the new formulation from [70], resampling occurs after the measurement 
update of both PF and KF. Therefore, the quantities computed for the mea- 
surement update of the PF can be reused for the KF measurement update. In 
particular, each particle is only evaluated once in line 1 of each MPF iteration 
instead of twice as with the previous formulation in [155]. Adopted from [78]. 
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Algorithm 3: Marginalized particle filter derived from [70, 155]. Adopted 


from [78]. 
. AL AN L- pL oN L N 
Input: Aj, Aj Cp-1, P p p Rp Qp Rp-1 Up Up, 
aL,- PpP „N,-,p 
*p-1 en Y p-1 
/* la PF measurement update x/ 
For p = 1,..., P, compute the particle importance weights de using the 


11 


12 


13 


14 


likelihood g? = N (SP, Sp), $” = CaP aha +c (ah 


SAS CaP Cl + R,-ı and compute the normalized weights 
P 


~P _ Ip 
Ip = TP 
p’=1 4P 
/* 1b KF measurement update x / 
2 L,+,p 2aL.-,p L-pT el ( - $P) 
a Ly +P, 162-18,-1 Yp-ı =) 
L+ ore ar Le 
EA = Por Pra lp a P p 
/* lc Resampling x / 


Resample P particles with replacement, 


probability (2749 = 274P) =g. 

fi — mean of 2572, p=l,...,P 

/* 2a KF time update x / 
ip ALR? + uk 

È, E AR + un 

PL- e ALPHY (AL) +at 


PE ANP ES (Ap) +07 
pél- Ape (ab) 
PiE- (PER) 


/* 2b PF time update x/ 
For p = 1,..., P, predict new particles, ne ~N Ge Pi). 
/* 2c Mixing step, update KF x / 


ahr? a shor wht (pp) (ajos Ez) 
Po — pL- pre- (P5) pee 


a EI pa NS 
Output: PE a P p 
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3.4.3 Nonlinear recursive B-spline approximation 
algorithm 


O Nonlinear recursive B-spline approximation (NRBA) adapts a B-spline func- 
tion f(s) with degree d iteratively to the data set from Section 3.2. Algorithm 4 
states the instructions for one iteration of NRBA, which is based on the MPF. 


In each iteration p, NRBA modifies f in J € N consecutive spline intervals. 


x : aL ¿E A . : 
Each linear particle or = Com 4 as al and each nonlinear particle 
„N, aN a A a y f 

X) P = es = a A contains estimates for J = d + I control points of 
ff. Kp = (Kp;> Kpa +++>Kp,) denotes the knot vector comprising K = J+d+1 


knots. The resulting definition range D, of f is given by Dp = [Kpa Kpy+1)- 
NRBA checks if sp is in the definition range of the previous time step, Dp-1. 
If not, D,-ı needs to be shifted such that sp € Dp. A shift can be conducted 
in the MPF time update. The result of the time update is the a priori estimate 
x7, In the following measurement update, sp is needed again to compute the 
measurement matrix Cp, and then, to take into account y y The result of the 
measurement update is the a posteriori estimate Lo 


Figure 3.5 depicts the allocation of available data points and computed estimates 
x to KF iterations in RBA versus MPF iterations in NRBA. The arrows indicate 
the needed information for computing the estimates. The KF is initialized with 
% and conducts in each iteration a time update first and then a measurement 
update. Therefore, P iterations are required for P data points. In contrast, 
the MPF performs the measurement update first and is initialized with £5. 
Therefore, y > has to be stored and sp, Sp+1, and y Pp have to be provided for 
iteration p+ 1. Hence, in order to take into account all data points, an additional 
iteration compared to the KF is needed. By definition, (s1, y,) is used for 
computing X) and s, for &,,,|, as indicated by the dashed arrows. 
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Algorithm 4: Nonlin. recurs. B-spline approximation. Adopted from [78]. 


ANP 2 L,- Bug 
Input: Kp- 1 A y PE Rp Sp Sp-1 Y p Y po Kp P 
IRA 
J < length ($7); K — length (xp-1); d- K-J-1; I- J-d 
/x Quantities for MPF measurement update x/ 


Compute u such that s,-ı € [Kp-1,,> Kp-1 
Cp-1 € RY"! from (3.47) 
/* Quantities for MPF time update x / 
a<-0 
if sp > Kp-1,,, then 
if sp > Kp-1, then 

| oa-d+l 
else 

| Compute o such that s, € 
end 
else if sp < kp-1,,, then 
if sp < Kp-1, then 

| ao -(d+ 1) 
else 

| Compute o such that sp € [Kp-1 a0 Kp-luo) 
end 


); Vp-1 — length (y ,_ı) 


H+1 


sen kolor) 


end 
ifo > 0 then 
Ke last element of tro 


Kp, Us, uN from (3. 18), 3 49) and (3.53) 
else 

Ke fst element of $i 

Kp, Up, un from (3. 18), 3. 55) and (3.56) 
end 


Compute u such that sp € [Kp Kp,,,,) 

H u+ 
AL, QE, AN and QY from (3.48), (3.51), (3.52) and (3.54) 
[PH na er gh 2] Algorithm 3 (AL AY ¡Cha Pa i; 


p >"p-V p 


L N R L $ aL,-,p „N,-,p 
Q7, Q Rp-1, Up» 2 Sp 
7 e+ aL,- =P pL,- 
Output: Kp, 5 8p ‚P 
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Initialization 
Each linear particle iP is initialized with a = zinitj jx], and each non- 
linear particle ar is initialized with 

„N,-,p _ -hit = 

Xo =X 1,x1 + chol (PI zxJ) -rmdjy1. 


Hereby, 17xı is a J x 1 matrix of ones and x!" indicates an initial value equal 
to the scalar measurement yy y referring to f. chol(-) computes the Cholesky 
factorization, and rnd ,xı is a J x 1 vector of random values drawn from the 
standard normal distribution. The covariance matrix of a priori estimation error 
of linear states, PP”, is initialized with Po = PI xy. [ y: denotes a J x J 
identity matrix. 


The large scalar value p causes X, to quickly change such that f adapts to the 
data. Provided that the values in Q% are small, the values in Pas decrease as p 
grows according to line 8 of Algorithm 3. Small elements in PE correspond 


E A : „L,-, 2N,=, 
to certain estimates. Therefore, the particles £, P and & p 2 


are slower to 
be updated using measurements such that f converges. Analogous statements 


hold for PET according to line 9 of Algorithm 3. 


Hence, the process noises are defined as Q% = GH ¡xj and Qay = NI jx 
with small positive values g/ and GN, respectively. 


Measurement update 
The measurement update from line 1 to line 4 of Algorithm 3 adapts f(s) 
based on (sp-1,Yp-1). The v-th dimension of y ,,_, refers to either f itself or 


a derivative of f. Therefore, the v-th row of the Vp- x J measurement matrix 
Cp-1 reads 


o” 
Cp-1ya, = (oran dad): (3.47) 
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Time step 1 2 3 n 
32% + 1— p 
Data points (s1,¥1) (s522) (53,3) (Sn, Yn) 
: 2= at a= at a at a- at a= a+ a- 
Estimates I HH HH HH Xn n Kar 
—y— ë — ra — y 
KF Init. Iter. 1 Iter. 2 Iter. 3 Iter. n 
A N rm oy 
MPF Init. Iter. 1 Iter. 2 Iter. 3 Iter. n + 1 


Figure 3.5: Allocation of availabe data points and computed estimates X to KF iterations in RBA 
versus MPF iterations in NRBA. Arrows indicate the needed information for computing 
the a priori estimates X” and the a posteriori estimates +. By definition the MPF uses 


(s1, y,) for computing x and sy for computing X, ,,. Adopted from [78]. 


whereby s,-ı € [Ky, Ku+1) and r € No. Algorithm 4 computes C ,-ı in line 3 
using (3.47). 


The value of the nonlinear measurement function c depends on the nonlinear 
particles > . Furthermore, c can depend on additional quantities that vary 
with the application and are not stated in Algorithm 3. 


The diagonal Vp x Vp covariance matrix of measurement noise R,_ı enables a 
relative weighting of the dimensions of y „_, because the influence of the vth 
: : A AL AN, 
dimension of the measurement error e}, = ( Yn-ı- IP ) on °°? and? 
P P p-1 p-1 


decreases with a growing positive value Rp-1,,.,. 


Time update with shift operation 
Based on a comparison between K p-1 and sp, NRBA decides if a shift operation 
of the B-spline function definition range is required to achieve that s, € Dp. 


The variable o calculated from line 4 to line 17 of Algorithm 4 states the shift 
„L,-,p 


direction of D,-ı and by how many positions components in Kp-1, 1 and 
2N-"P need to be moved for that purpose. o > O indicates a right shift of Dp-1, 


p-1 
o < 0 indicates a left shift, and o = O means that no shift is conducted because 


Sp € Dp-1- 
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Algorithm 4 expects that, for o > 0, the |o| additionally needed knots are the 
o last entries of the knot vector Kp = (Kp;» Kpa» - - -» Kpg ) and that they are the 
—o first entries of Kp fo < 0. 


Case 1: Right shift of definition range (c > 0) 


The new knot vector xp is given by (3.18) and line 6 of Algorithm 3 updates 


a to ER using the state transition matrix 
p 
Al, = Ap from (3.19) (3.48) 
and the input signal vector 
u), = u, from (3.21). (3.49) 
„L,+,p 


Thereby all entries of x, are moved to the left and the last o entries of 
AoE have an arbitrary initial value x: 


: 
BE Hg) (3.50) 


or’ J- 


„L-,p _ (25 aL 
Xp Moto *p-1 


During a right shift of the definition range, x is set to the last element of £59» 
which is determined during the preceding call of Algorithm 3 in line 5. This is 
based on the assumption that a) is a good initial value in the magnitude of 
the data. 


Additionally, line 8 of Algorithm 3 updates Po to Ppr using (3.48) and 


B ifh=g^Q 
Q; ER’ with Q., =43%, ifħ=g^-0 (3.51) 


0, otherwise 


. . . . L,+ 
with O given by (3.24). The update operation moves the elements in Ppa to 


the top left and replaces the zeros on the last o main diagonal elements of a 
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with p in order to get large values on the last o main diagonal elements of P 
and a fast adaption of the initial estimates x to the data points. 


In line 7 and line 9, Algorithm 3 computes the quantities Es and PET that are 
needed for the PF time update. The calculations of the state transition matrix 
AN with 

AN = A, from (3.19) (3.52) 


and the input signal vector uN with 


un = u, from (3.21) (3.53) 


are analogous to those for the linear quantities. Q” uses gy instead of gr: 


B ifh=grQ 
QY ER’ withQN „= 17", ifh=g1-0 (3.54) 


Peg;h 


0, otherwise. 


Case 2: Left shift of definition range (c < 0) 


The updated knot vector is given by (3.25), the input signal vector for linear 
states u? reads 
ul, = u, from (3.26) (3.55) 


and the input signal vector for nonlinear states u is given by 


un = u, from (3.26). (3.56) 


Additionally, x is set to the first component of a ie 


Since AL and AY are identical in the chosen state-space model, computational 
effort can be saved when calculating the covariances and cross-covariances from 
line 8 to line 11 in Algorithm 3. Adopted from [78]. 
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3.4.4 Effectiveness of nonlinear recursive B-spline 
approximation 


O Numerical experiments conducted in [78] demonstrate the effectiveness of 
Algorithm 4 in comparison with solutions determined using the LM [115] al- 
gorithm. Therein effects of the number of simultaneously adaptable spline 
intervals and the particle count on the NRBA solution are also investigated. A 
MATLAB implementation of the experiments is provided in [76]. 


General experimental setup 


The data set ((sp, Y») p=1,2....P 1s defined according to Section 3.2, whereby 


Sp =0.25+0.5-(p=1, 


40, if 80 < sp < 120 
y = 
ý 30, otherwise (3.57) 


Yp,2 = Yp,3 = Yp,4 = 0 Yp 
P = 400. 


A B-spline function f(s) with knot vector k = (-30,-—20,...,230) and de- 
gree d = 3 approximates the data. Thereby, it is supposed that y,,ı refers 
to f, yp,2 to the first derivative of f, yp,3 to the second derivative of f, 
and yp,4 to the value of the nonlinear measurement function c, which is 
defined as a quadratic B-spline function with k = (-5,0,...,70) and x = 
(0, 0, 0, 0.25, 1.5, 5, 5, 0, 0, 6, 8, 8,8)”. 


c depends on the value of the approximating function f(s) and is displayed 
in Figure 3.6. The input variable f(s) of c is restricted to the definition range 
[5, 60] of c. 


The diagonal measurement covariance matrix Rp € RA with Rpa = 1, 


Rp, =5: 10, Rp = 5° 107° and Rp44 = 0.8 or 10°, respectively, comprises 
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c (f(s) 


oN KD OO 


Fs) 
10 20 30 40 50 60 


Figure 3.6: Nonlinear measurement function c (f(s)) depending on the value of the B-spline 
function f (s) approximating the data. Adopted from [78]. 


the reciprocal weights of y,,1, Yp,2, Yp,3 and yp,4. The reciprocal weight values 
for the first three dimensions of y,, avoid that f oscillates and cause that f 
smooths the jumps in the first dimension of the measurements. With Rp44 = 
0.8, the nonlinear target criterion c (f(s)) = 0 is weighted strongly, whereas it 
is almost completely neglected with R p44 = 10°. 


Depending on the applied algorithm, solutions for the former weighting case 
are denoted by NRBAN or LMN indicating the nonlinear problem. Conversely, 
solutions for the latter case of a quasi-linear approximation problem are denoted 
by NRBA! or LM". 


Solutions for two different numbers of spline intervals J are analysed. For J = 1, 
K is initialized with Kọ = (—30, 20,...,40) which leads to an initial definition 
range [0, 10) of f(s). For I = 3, x is initialized with ko = (—30, 20,..., 60) 
and the resulting definition range is [0, 30). In both cases NRBA approximates 
the data by repeatedly shifting the function definition range to the right. Each 
time, an additional knot value k,,, needs to be provided in the vector Kp. For 
I = 1 these values are kp, = 50,60,...,230 and for J = 3 they read kp, = 
70, 80, ...,230. In order to display the NRBA results for the whole data set, all 
values that are moved out of NRBA matrices and vectors are stored elsewhere. 


The remaining NRBA parameters are 7" = 0.005, GN = 0.25 and p = 30. The 


LM algorithm uses x” = 30 as the initial value for each control point. 


Due to the included PF, NRBA is a nondeterministic method. For each setting, 
50 runs are performed and for each run the normalized root mean square error 
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(NRMSE) between the B-spline function determined by NRBA, /nrBa; and 
the B-spline function according to LM, fim, is calculated as follows: 


1 


maxp=1,...P{fLM(Sp)} — minp=1,..,P{fLmM(Sp)} 


\ EP, (fsal) - fum(sp)) 
P 


NRMSE = 


(3.58) 


The terms NRBA™*¢ and NRBA" refer to the NRBA solution with the median 
or maximum NRMSE, respectively, in each set of 50 runs. 


Results 


Figure 3.7 shows for both the quasi-linear (L) and the nonlinear (N) problem 
the approximating functions NRBA”“ and NRBA™* compared to the LM 
solutions. Black dots depict the first component y, of 40 of the 400 data 
points (sp, Yp). Dashed vertical lines indicate knots, whereby the first and last 
knots are not shown. 


Figure 3.7a shows NRBA approximations for / = 1 and P = 6561 = 94. In this 
case the MPF state vector comprises four linear and four nonlinear components 
and the PF creates nine samples per nonlinear state dimension. 


At f(s) = 30, the deviation between the value of c and its target value y, 4 = 0 
has a local maximum (c.f. Figure 3.6). The nonlinear problem penalizes this 
deviation strongly; hence, NRBAN and LMN avoid f(s) = 30. In contrast, 
NRBA! and LM! approximate data with Yp,1 = 30 closely. 


Data and knot vector are symmetrical to the straight line given by s = 100. 
Since the LM algorithm processes all data simultaneously in each iteration, the 
solutions LM" and LMN in Figure 3.7a reflect this symmetry. 


In contrast, NRBA processes the data from left to right and can only adapt 
some control points at a time. For J = 1, these are the four control points that 
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Figure 3.7: B-spline functions f with NRBA for interval counts J and particle counts P and with 
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linear (L) and nonlinear (N) problems. Only 40 data points (sp, yp,ı) and a subset of 


median (med) and maximum (max) NRMSE compared to the LM solution for quasi- 
knots x is shown. Adopted from [78]. 


3 Data approximation with B-spline functions 


influence the B-spline function in the spline interval in which the current data 
point lies and for / = 3 additionally the two control points that affect the two 
spline intervals to the left. 


NRBA! and NRBAN are both asymmetrical and mostly delayed with respect 
to LM! and LMN. For NRBAN, the asymmetry is less distinct because the PF 
removes states that translate to a large delay more quickly from the particle set 
because they create a larger error. Additionally, the range of values in NRBAN 
is smaller than in NRBAT so that a present lag is less obvious. 


For the same weighting NRBA™¢ and NRBA™* differ only slightly, which 
suggests that, for the given setup, P = 6561 suffices for a convergence of 
NRBA solutions. 


With RBA for linear weighted least squares approximation similar numerical 
experiments but without any nonlinear approximation criterion were performed 
in Subsection 3.3.4 and [80]. For J = 1, a strong asymmetry and delay are 
observed with RBA, analogous to NRBAT in Figure 3.7a. The delay decreased 
as I was increased because the filter was able to update more control point 
estimates with hindsight based on PF*. 


By increasing J to three with NRBA, the dimension of the state space also 
increases to six linear and six nonlinear components. The PF in NRBA samples 
the state space less densely, unless the particle count is increased exponentially 
with 7. 


Figure 3.7b displays the results for the quasi-linear approximation problem for 
the case of keeping the sampling density per nonlinear state space dimension 
constant by choosing P = 625 = 5* for I = 1 and P = 15625 = 5° for 
I = 3. Overall the NRBA solution is more symmetrical with J = 3 but a 
comparison of NRBA”“ for J = 1 and J = 3 indicates that the delay for 
s > 120 is not reduced. Instead, the ability to adapt more control point estimates 
simultaneously sometimes leads to undesirable results, e.g. the too low course 
between s = 40 and s = 60 as well as the overcompensation of the delay 
between s = 60 and s = 75. 
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For J = 1, NRBA™* differs more from NRBA”“ and shows larger oscillation 
amplitudes than for J = 3. This suggests that P = 625 is not sufficient for 
a convergence of NRBA for J = 1. However, even with 625 particles for 
I = 1, the required increase to P = 15625 for J = 3 is already quite strong. 
Keeping the sampling density constant quickly becomes infeasible, especially 
if computation time constraints are present [70]. 


All other factors held constant Figure 3.7c shows the results for the nonlinear 
approximation problem, which support the previously drawn conclusions. Ad- 
ditionally, the conflicting target criteria in the nonlinear approximation problem 
cause a larger stabilization period at s < 20. 


Figure 3.7d depicts the effect of increasing / from one to three while main- 
taining the particle count of Figure 3.7a. For J = 3 NRBA™* differ much 
more from the corresponding N RBA". This indicates that more particles are 
needed for convergence for J = 3. Furthermore, for J = 3 these differences are 
much larger for NRBAN than for NRBAL. 


A comparison between corresponding NRBA™ solutions in both figures shows 
only a small approximation improvement from increasing 7 for the chosen 
setup. In Figure 3.7d NRBAN temporarily decreases below f(s) = 30, the 
position of the maximum of c (c.f. Figure 3.6). This is a case in which NRBA 
chooses a locally but not globally good solution. More detailed investigations 
on the effects of particle count P on convergence can be found in the cited 
prepublication. Adopted from [78]. 


3.5 Scientific contribution 


This chapter presented two novel algorithms for the weighted least squares 
(WLS) approximation of a set of data points with a B-spline function. These are 
the algorithm recursive B-spline approximation (RBA) for the case of a linear 
WLS problem and the algorithm nonlinear recursive B-spline approximation 
(NRBA) for the case of a nonlinear WES problem. 
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Unconstrained WLS B-spline function approximation 


Linear 


Nonlinear 


ms] [65%] a] [iw] em] 


Figure 3.8: Classification of unconstrained weighted least squares (WLS) B-spline approximation 
problems and suitable algorithms. 


Table 3.1: Comparison of B-spline approximation methods. *: Only applies to WLS. **: Also 
applies to LM if the sliding window method in [38] is used with a window size very 
small compared to P. Adopted from [80]. 


LS/LM WLS/LM 


Feature Cue call) (multiple calls) RLS/KF RBA/NRBA 
Number of sabl 

data points er 5 Bounded Unbounded Unbounded Unbounded 
Time complexity O(P)* O(P)** O(P) O(P) 
Approximation interval Fixed Variable Fixed Variable 
Determination of total 

number of control points At. During At During 
being estimated beginning run-time beginning run-time 


Figure 3.8 subdivides the unconstrained WLS B-spline approximation problem 
and illustrates the problem types to which RBA, NRBA and the well-known 
algorithms WLS, RLS, KF and LM can be applied. 


Table 3.1 compares features of the mentioned approximation methods. With 
LM the computational effort in each iteration depends on the number of data 
points P. However, the number of iterations performed depends on the specified 
tolerance and cannot be deduced solely from P. Therefore the time complexity 
statement in Table 3.1 does not apply to LM in the general case. However, the 
sliding window LM approach in [38] has time compexity O(P) if the window 
size is very small compared to P. 
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3.5 Scientific contribution 


The design of RBA, NRBA in this chapter as well as their applications in the 
following chapters are limited to the B-spline representation and functions. 
However, the limitation to functions and B-splines in this work did not lead to 
a loss of generality and the approach of RBA and NRBA can be generalized 
to curves and surfaces of all spline representations that offer local control, so 
that the benefits of RBA and NRBA stated in the following can be directly 
transferred to them as well. 


This is because according to Section 2.3, curves and surfaces are direct exten- 
sions of functions. Moreover, spline representations with local control only 
differ in the geometry matrix M and therefore in the blending functions, but 
they all share the same structure of (2.17) for curves and of (2.19) for surfaces. 
Their local control results from the bounded interval, in which a blending func- 
tion is nonzero. As a result, other local spline representations suffer from the 
same identified research gap. 


Recursive B-spline approximation 


O The RBA algorithm prepublished in [80] solves a linear WLS approximation 
problem iteratively using a KF which estimates the control points of the B- 
spline function sequentially. Therefore the total computational effort increases 
linearly with the number of approximated data points. 


The main contribution is to use the time update of the KF for a shift of estimated 
B-spline control points in the KF state vector in combination with a shift in the 
B-spline knot vector. The shift operation enables to shift the definition range 
such that it is always possible to take into account the latest data point for the 
approximation. 


Thereby RBA overcomes the limitation of other recursive algorithms based on 
the KF that can only approximate data points within the initially chosen and 
fixed approximation interval. 
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RBA is especially advantageous in online applications in which the magnitude 
of the data points is not exactly known or changes over time because then data 
points outside the initially chosen bounded B-spline function definition range 
can occur. 


RBA is also beneficial when a tradeoff between low computational effort and 
high approximation quality is needed because the shift operation of RBA allows 
to reduce the size of the KF state vector in both online and offline applications. 
A smaller state vector causes less computational effort. 


Numerical experiments in [80] and Subsection 3.3.4 show that the RBA re- 
sult converges to the WLS solution as the size of the state vector is increased. 
Additionally the experimental results reveal that few simultaneously adaptable 
spline intervals suffice for good approximation results. 


The number of required shift operations can increase when the size of the 
state vector is reduced to lower the computational effort. The experiments 
indicate that shift operations influence the control point values of the resulting 
approximation function only in a magnitude that is close to numerical accuracy. 
Each shift operation comes at the expense that a part of the approximation 
result is forgotten in order to keep the sizes of matrices and vectors constant. 
A growing approximation interval can be realized by storing matrix and vector 
elements separately from RBA before they are overwritten. Adopted from 
[80]. 


Nonlinear recursive B-spline approximation 


O The NRBA algorithm prepublished in [78] is a generalization of RBA for 
NWLS problems which result from target criteria that depend on the control 
points nonlinearly. 


NRBA determines a B-spline function such that it approximates an unbounded 
number of data points with respect to both linear and nonlinear target criteria. 
The approach uses a MPF for solving the approximation problem iteratively. 
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In the MPF, a PF takes into account target criteria that do not relate to the 
control points in a linear fashion whereas a KF solves any linear subproblem 
optimally for each particle [70]. As the values of the B-spline function and its 
derivatives depend linearly on the control point values, linear target criteria will 
occur in most approximation applications. 


The MPF can take into account the exactly known values of the B-spline basis 
functions and does not need to estimate them like most other nonlinear filters 
do. Taking advantage of the linear substructure of the problem allows to reduce 
computational effort and achieve better results compared to purely nonlinear 
filters like a PF [155]. 


The features and benefits of RBA that result from the shift operation also apply 
to NRBA. Numerical experiments in [78] and Subsection 3.4.4 investigated 
the effectiveness of the approach in comparison to the LM algorithm and il- 
lustrated the effects of selected NRBA parameter values on the approximation 
result. Provided that the parameters are chosen appropriately, the solution of the 
proposed method is close to the LM solution apart from a slight filter-typical 
delay. 


NRBA use cases are NWLS problems in which a linearization of nonlinear cri- 
teria is not desired or promising, for example because of distinct nonlinearities. 


For linear WLS problems RBA should be used instead of NRBA. RBA is based 
on the KF, which computes an optimal solution [205]. For linear problems 
NRBA can at best reach the same approximation quality provided that the 
particle count is large enough, which requires more computational effort. 


Furthermore, with NRBA the approximation depends more strongly on the 
parameterization of the underlying filter algorithm than with RBA. 


Increasing the number of control points that NRBA can adapt simultaneously 
is not as unambiguously beneficial for the approximation result as with RBA. 
Moreover, with NRBA the increase of control point count usually needs to be 
combined with an exponential increase of the particle count in the PF for an 
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Figure 3.9: Exemplary use cases for different operating modes of RBA and NRBA. Optimization 
range (OR). 


improvement of the approximation. This is known as curse of dimensionality 
and a general problem of sampling-based nonlinear filters. Adopted from [78]. 


Use cases in terms of operating mode 


Use cases of RBA and NRBA can be classified according to whether the knot 
vector changes only with a shift operation or with each iteration: 


e Knots change only with each shift operation: This mode assumes that the 
independent variable will not take values large enough to create an arith- 
metic overflow during operation. Therefore the knot vector changes only 
with a shift operation. Figure 3.9a illustrates this mode at the example of 
trajectory optimization, whereby the independent variable 7 refers to the 
planned trajectory time. From the top plot to the bottom plot the compu- 
tation time progresses and the number of processed data points increases. 
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In the example, ¢ grows with the data, but a limited planning horizon, 
that does not lead to arithmetic overflows, can be assumed. Ideally, t 
increases much faster than real time so that there is no need to evaluate 
the trajectory in the optimization range (OR), in which it can change. 
Chapter 5 will investigate this use case in more detail. The examples 
in Subsection 3.3.4 and Subsection 3.4.4 also belong to this operating 
mode. 


Knots change with each iteration: This operating mode is required in 
case of an unbounded independent variable, that eventually will take a 
value that causes an arithmetic overflow. Figure 3.9b illustrates this mode 
at the example of an analytic representation of a signal or its derivative or 
its integral over the recent past up to the present point in time. From the 
top plot to the bottom plot a new measurement is received with each plot. 
If the independent variable t were assigned to the observation number, 
eventually an overflow in the knot vector would be encountered. Let 
dt denote the time interval between the two latest measurements. For 
simplicity, let dt be constant. An overflow can be avoided by subtracting 
dt from the knot vector after each iteration. Thereby the latest measure- 
ment is always assigned to ¢ = O and the knots are continuously shifted 
to the left in Figure 3.9b as their values keep decreasing until they are 
removed from the knot vector and the OR of RBA or NRBA. The signal 
representation application requires to evaluate the function within the OR 
unless a time delay of the magnitude of the OR is acceptable and a some 
knots and control points can be saved when they leave the OR. If so, the 
function can be evaluated outside of the OR and further into the past. 


Differentation from previous works 


In the previous project e-generation [22, pp. 29-35] an iterative data approxima- 


tion by a functional representation was developed as well and was patented by 


F. Bleimund and S. Rhode [21]. This previous approach has in common with 


RBA that it also solves the linear WLS optimization problem (3.13) iteratively, 
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whereby the target criteria refer to the value of the functional representation 
and its first two derivatives. Moreover, the optimization problem is solved 
iteratively using a KF and a shift matrix. 


However, in the previous project, a polynomial as in (2.8) is adapted to the data. 
This polynomial exhibits the features discussed in Section 2.3 along with its 
flaws compared to a spline representation, such as coupling between degrees 
of freedom and polynomial degree, risk of oscillations from the Runge’s phe- 
nomenon, global effect of coefficient value changes and geometrically difficult 
to interpret coefficients. 


Due to the global control of the polynomial, a shift operation and therefore the 
KF would not be needed if the approximation were performed analogously to 
the case depicted by Figure 3.9a, in which the value of the independent variable 
can take large values. For this case, RLS would suffice. 


However, the approach assigns the latest oberserved data point always to the 
same value of the independent variable, analogous to Figure 3.9b. This is 
achieved by defining the shift operation via an upper triangular matrix such that 
a coefficient vector is calculated for a polynomial that is shifted by dt to the left 
with respect to the current polynomial during the time update of the KF. This 
shift operation for the polynomial corresponds to the subtraction of dt from the 
spline knot vector as described above for the second operating mode. 


In contrast, the shifting operation stated in the previous sections of this work is 
needed in general for both of the above described operating modes. It modifies 
both the knot vector and the control point vector and occurs depending on the 
knot vector and the independent variable of the latest data point. 


The algorithms presented in this work are superior to the one developed in 
e-generation in that they allow to benefit from the advantages of spline repre- 
sentations with local control during an iterative data approximation. The NWLS 
problem adressed by NRBA was not considered in the previous project. 
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4 Models of research vehicle and 
reference route 


This chapter describes models of the research vehicle and the reference route. 
The content of this chapter and the previous one will then be combined for 
trajectory optimization in Chapter 5 and ALC in Chapter 6. 


The research vehicle is an all-wheel-driven BEV based on the Porsche Boxster 
(type 981) which was developed during the research project e-generation. Fur- 
ther information about the vehicle is provided in the remainder of this chapter 
as well as in [14, 16, 207]. 


Detailed nonlinear models allow a representation of the dynamic vehicle behav- 
ior that is close to reality. Such models take into account the kinematics of the 
vehicle and of its subsystems like the wheel suspension. They also comprise 
models of tire forces. However, a model should also focus on the aspects that 
are important for the actual application and simplify where possible. Depend- 
ing on the desired level of detail and the application, different vehicle models 
are known in literature [157, pp. 5-6]. 


The purpose of the vehicle models in this work is to investigate the energy 
consumption of the research vehicle under ALC on a reference route. According 
to [187, p. 28] with a normally experienced driver the absolute value of the 
lateral acceleration is below 3.5 m/s? and the longitudinal acceleration ranges 
roughly from —2.4 m/s? to 1.8 m/s”. For comfort and acceptance reasons, an 
ALC should not exceed these limits most of the time, therefore models for 
low dynamics suffice. Simplifying assumptions of this work for longitudinal 
dynamics include: 
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Figure 4.1: Lateral view of research vehicle at slope of angle œ with acting climbing force Fa, 
rolling resistance Fon, cornering resistance Foor, air resistance Fy;, and inertial force 
Finert- Vehicle drawing from [46]. 


Lifting, rolling and pitching motions of the vehicle are neglected. 
e The vehicle mass is concentrated in the center of gravity (COG). 
+ The vertical wheel force, also known as wheel load, is constant. 
e Longitudinal tire slip is neglected. 


For describing the lateral dynamics the linear single track model is used. 


Section 4.1 states the coordinate system. Section 4.2 describes relevant driving 
resistances and Section 4.3 the power train of the research vehicle. Section 4.4 
explains how the resulting energy consumption for a given route can be de- 
rived and reasons the approach for optimization of energy consumption that 
the ALC takes. Section 4.5 describes vehicle models for ALC tests in a simu- 
lation environment and for investigations of the resulting energy consumption 
of the research vehicle. The reference route is mentioned in Section 4.6 and 
Section 4.7 derives two simplified adaptive vehicle models for use within the 
ALC. 
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4.1 Coordinate system 


4.1 Coordinate system 


Several coordinate systems are used to describe vehicle position and orientation 
[157, pp. 17-31]. This work uses a body-fixed right-hand vehicle coordinate 
system with origin in the vehicle COG and axes x, y and z. The axes are 
perpendicular to each other and shown in Figure 4.1 and Figure 4.2. The x axis 
is the vehicle longitudinal axis and points towards the vehicle front, whereas 
the y axis is the vehicle lateral axis that points towards the left vehicle side. 
Since rolling and pitching movements of the vehicle with respect to the road 
surface are neglected, the plane that x and y axes span is parallel to the road 
surface. The z axis is the vertical axis that points upwards and is perpendicular 
to the road surface. 


4.2 Driving resistances 


This section states the main forces acting on a vehicle and the force equilibrium 
that they form at the wheels. Figure 4.1 and Figure 4.2 illustrate these forces. 


4.2.1 Climbing force 


O The climbing force Fa results from the gravitational force acting on the COG 
and is given by 
Fa = Mynel + 8: sin(a). (4.1) 


Fy depends on the vehicle mass myhe and road slope angle œ which is measured 
between the road surface and the plane that is perpendicular to the direction of 
gravitational force myncı - g With gravitational constant g. Adopted from [79]. 
The road slope angle œ can be computed from the road slope y stated in percent 
using 

a = arctan(y/100) [119, p. 95]. (4.2) 
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4.2.2 Wheel resistance 


The wheel resistance consists of several components, of which only rolling 
resistance and cornering resistance will be considered. Resistance components 
neglected in this work include the toe-in resistance and flood resistance. The toe- 
in resistance is caused by not perfectly parallel wheels and the flood resistance 
occurs on wet roads when the tire evacuates water [119, pp. 11-19]. Moreover, 
the longitudinal tire slip resistance [165] and the ventilation resistance of the 
rotating wheel [182] are not considered. 


Bearing friction as well as residual braking torque from the hydraulic brakes 
can also be assigned to wheel resistance but will be considered as resistances 
of the power train in this work. 


Rolling resistance 


O Rolling resistance results from damping forces of the deformed tire rubber. 
When driving straight on a dry road, the wheel resistance almost entirely comes 
from the rolling resistance Fon [119, p. 11] with 


Fron = Fr Myhel + 8 + cos(a). (4.3) 


Adopted from [79]. The normal force Myhe - g - cos(a) equals the sum 
of vertical front axle wheel force Fw: EA and vertical rear axle wheel force 
Fwnzra- In the static case the distribution of the vertical forces between front 
axle (FA) and rear axle (RA) can be determined from the wheel base (/), the 
distance between FA and COG (Ira) and the distance between RA and COG 
(Ira). The road slope influences the static load distribution between the axles 
as well but since the distance between COG of the research vehicle and road 
surface is small, this effect of road slope is neglected. 
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4.2 Driving resistances 


The rolling resistance coefficient fy depends on the vertical wheel force, tire 
pressure and slightly increases with velocity [24, pp. 50, 52-53]. For most 
passenger vehicles f; can be assumed between 0.007 and 0.014 [10]. 


Cornering resistance 
In curves the centrifugal force Foen with 


Feentr = Myhel * VVhcl,y (4.4) 


acts on the COG. The lateral acceleration Ýynci,y is given by 
. = 
Wohely = Vine * K (4.5) 


for small side slip angles and constant velocity [157, p. 228]. The road curvature 
K is the reciprocal of the curve radius. 


The single track vehicle model represents the fundamental driving dynamics 
for lateral acceleration smaller than vvyhe1 y = 4 m/ s? on dry roads. The model 
1s very idealized and uses the following simplifications: 


Lifting, rolling and pitching motions of the vehicle are neglected. 


The vehicle mass is concentrated in the COG. 


The vehicle velocity is assumed to be roughly constant, 1.e. quasi- 
stationary. 


Wheels on the same axle are represented by a single wheel in the center 
of the corresponding axle. 


The vertical wheel force is constant. 


Pneumatic trail and aligning torque resulting from the slip angle of the 
tire are neglected. 


Longitudinal tire forces are neglected [157, pp. 225-226]. 
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Figure 4.2: Top view of research vehicle with acting air forces. Environmental wind with velocity 
Vwind and angle Twina as well as vehicle velocity vynci result in relative air flow with 
velocity Vre] and angle 7;e¡. This causes longitudinal air resistance Fyjr, and lateral air 
resistance Fair y acting on the pressure point (PP). COG denotes the center of gravity. 
Vehicle drawing from [46]. 


In order to stay on the road, the tires need to provide lateral forces that counter- 
act the centrifugal force. These lateral forces are distributed among the axles 
depending on the position of the COG and created using tire slip angles. 


These slip angles cause the cornering resistance Foo, acting in longitudinal 
direction. Under the assumption of small road curvature, small steering angle, 
small slip angles and identical tires, Foo, according to the single track model 
reads 


Beyi 1 2. py 2 
Foor = A jE, Myhel} + Ra, Mvyhcl (4.6) 
2-cC l l 


whereby c denotes the tire cornering stiffness [64, p. 147]. 
4.2.3 Air resistance 
The air resistance Fhir acts on the pressure point (PP) of the vehicle and results 


from the turbulences of air stream and friction of the air flow. The relative air 
velocity vre] and relative air flow angle Tye] both depend on the vehicle velocity 
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4.2 Driving resistances 


Vvnel, the environmental wind velocity vwina and the environmental wind angle 
Twind aS Figure 4.2 illustrates: 


= 2 2 
Vrel = fig + Vivina + 2 + Vvhel * VWind * COS(Twind) 
(4.7) 


; VWind . 
Tre] = arcsin | —— | - sin(Twind) 


Vrel 


The angles are measured anticlockwise starting from negative direction of the 


x axis. O The longitudinal air resistance Fair x is given by 


Fatx = Cx (Tra) AS Vee (4.8) 
Cx is the dimensionless longitudinal drag coefficient that describes the shape 
of the vehicle and depends on the relative air flow angle Tye. A denotes the 
effective vehicle cross-sectional area perpendicular to the vehicle longitudinal 
axis and p is the air density. Adopted from [79]. Analogous to the centrifu- 
gal force, the lateral air resistance Fairy needs to be balanced by lateral tire 
forces. The required slip angles in turn cause a driving resistance component 
in longitudinal direction. 


Neglecting environmental wind leads to vie] = Vvhel» Trel = O and Fairy = 0. 
Cx(Treı = 0) is also denoted cy [157, pp. 212-214]. For passenger vehicles 
Cw usually lies between 0.20 and 0.40 and A ranges from 1.5 m? to 2.5 m?. 
With increasing vehicle velocity Fair x becomes the largest driving resistance 
component [24, pp. 50-51]. 


4.2.4 Inertial force 


Inertial forces arise when vehicle velocity changes. Analogously variations 
of the angular velocity of power train components cause inertial torque. Each 
axle of the research vehicle is driven independently as Figure 4.3 illustrates. 
Due to the fixed gear ratio of each gear box and the assumption of negligible 
longitudinal tire slip, angular velocities are proportional to the vehicle velocity. 
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Computing the equivalent mass of the power train allows to summarize both 
inertial force and inertial torque in the inertial force Finet with respect to the 
vehicle longitudinal acceleration Vypci,x [62, pp. 16-17]: 


Finert = (Myncı + Meq,FA + Meq,RA) * Vvhel.x (4.9) 


Meg, F A 18 the equivalent mass of the power train components that drive the front 
axle FA and meq,ra denotes the equivalent mass of components that drive the 
rear axle RA: 


2 
JEM,FA ‘ig + Jog + 2+ Jwneel,FA + 2 * JBrake,FA 


Meg,FA 


r2 
dyn, FA (4 10) 


2 
JEM,RA ‘ig + Jog + 2: Jwheel,RA + 2: JBrake,RA 


Meg,RA 2 
dyn,RA 


Jem is the mass moment of inertia of the corresponding electric motor (EM) 
that is connected via a gear box (GB) with integrated differential to the wheels. 
Both gear boxes have the same gear ratio ig = 9.59 [14]. Jag is the mass 
moment of inertia of a GB, Jwheeı the mass moment of inertia of a wheel, JBrake 
the mass moment of inertia of a brake and rayn the dynamic wheel radius. 


4.2.5 Force equilibrium 


O The sum of the driving resistances from (4.1), (4.3), (4.6), (4.8) and (4.9) 
equals the traction force Fira. between all tires and the road surface: 


Frac = Fo + Fion + Fair t Feor + Finert (4.11) 
The mechanical traction power Pirac.mech Corresponding to Firac reads 


P trac,mech = Firac * VVhel- (4.12) 
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Figure 4.3: Model of power train of research vehicle. 


Providing Frac requires longitudinal tire slip. Under the assumption of negligi- 
ble slip losses, Pirac,mech equals the wheel power Pwheel provided by the power 
train [157, pp. 152-153]. Adopted from [79]. 


4.3 Power train 


Figure 4.3 depicts the power train of the research vehicle. In various compo- 
nents of the power train power losses occur. These power losses are measured 
during component and vehicle tests on test benches and are stored in look-up ta- 
bles, which have been provided for the research vehicle. Depending on the kind 
of input power, components of the power train can be divided into mechanical 
components and electrical components. 
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4.3.1 Mechanical components 


The wheel power Pwneeı provided by the power train reads 


Pwheel = TFA * @wnl,FA + TRA * @WhLRA (4.13) 


with front axle torque Tra, rear axle torque Tra, front wheel angular velocity 
@wni,ra and rear wheel angular velocity WwhLRA- 


The torque distribution between Tra and Tra is computed by the motor elec- 
tronic control unit (ECU) as a function of Frac. At each axle a gear box (GB) 
with open differential distributes the torque equally among the left and right 
wheel. 


When driving straight, wwniFa and wwnı,rA are given by 


_ VVhel _ VVhel 
@WhLFA = > @Whl,RA = 
Fdyn,FA Tdyn,RA 


(4.14) 


with dynamic front wheel radius rayn,ra and dynamic rear wheel radius rayn,RA- 
Between wheel and gear box undesired friction occurs in the wheel bearing, 
which results in braking torque. There is also residual brake torque caused by 
permanent residual friction between brake pad and brake disc. The model takes 
into account the front axle brake torque Tprx.Fa and the rear axle brake torque 
Terka. The values are derived by interpolating linearly with respect to the 
wheel angular velocities in look-up tables. 


At each axle a 1-speed gear box with gear ratio ig increases the motor angu- 
lar velocity and reduces the motor torque. Front axle motor angular velocity 
WeEM.FA and rear axle motor angular velocity WemraA are given by 


WEM,FA = 1iG*WWhLFA, WEMRA = İG * WWhLRA- (4.15) 
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Front axle motor torque 7km,ra and rear axle motor torque Tem. ra read 


_ Tra 
TEM,FA = Fra TBrk,GB,FA (TEM,FA, WEM,FA) » 
G (4.16) 
_ Tra 
TEMRA = TBrk.GB,RA (TEM,RA, WEM,RA) > 
G 


whereby TBrx GB.FA and Tgrk,ag Ra denote the braking torques occurring in the 
gear boxes during the not lossless power conversion. The absolute values of 
these quantities increase with the absolute values of motor torque and angular 
velocity and are determined using two-dimensional linear interpolation in look- 
up tables. 


4.3.2 Electrical components 


A permanently excited synchronous motor (PSM) with 120 kW drives the front 
axle and an asynchronous motor (ASM) with 144 kW drives the rear axle. The 
PSM has a higher efficiency whereas the ASM has lower drag losses. Therefore, 
mainly the PSM propels the vehicle and is supported by the ASM only in case 
of high power demands [14]. The motors are connected to the high voltage 
(HV) battery via power electronics (PE). The losses of both motors including 
the losses in the power electronics are modelled by two-dimensional linear 
interpolation with respect to motor torque Tem and motor angular velocity 
wgm in look-up tables. Ohmic resistance and magnetic resistance as well as 
switching losses contribute to the losses of EM and PE. Purac elec 18 the electrical 
traction power that results from the mechanical traction power Pirac,mech and 
includes the power losses between power electronics and wheels. 


While braking, the motors act as generators and convert mechanical traction 
power Pirac,mech into electrical traction power Pirac elec that is recuperated into 
the HV battery. For driving stability reasons, the majority of the recuperation 
power comes from the front wheels, which are connected to the more efficient 
PSM. When the braking power demand exceeds the recuperation capabilites, 
hydraulic wheel brakes are activated in order to fulfill the demand. Hydraulic 
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wheel brakes convert mechanical traction power into thermal power that cannot 
be reused for vehicle propulsion anymore. Additional hydraulic braking power 
is usually needed when braking at high velocity or when high deceleration is 
requested. Furthermore, during a standstill and at very low speeds hydraulic 
brakes are solely used. 


Due to the powerful drive train, the vehicle has a high recuperation power limit. 
Furthermore, simulations of the ALC for determining the energy consumption 
will not include standstills, high dynamic driving maneuvers and driving at high 
speeds. Therefore this work assumes that the vehicle can brake solely using 
recuperation, also denoted perfect recuperation by [62, p. 16], and does not 
consider hydraulic braking power in the vehicle model. 


Apart from the power electronics, the motors and the HV battery, a positive 
temperature coefficient (PTC) battery heating, an air conditioning (AC) com- 
pressor and a DC/DC converter for the 12 V circuit belong to the HV system. 
The on-board network (OBN) is fed via the 12 V battery or directly from the 
DC/DC converter. The on-board network (OBN) also powers a thermoelec- 
tric heat pump for heating the passenger cabin. Arrows in Figure 4.3 indicate 
possible power flow directions between electrical components. 


As a result the PTC power consumption Pprc, the AC power consumption Pac 
and the DC/DC converter power consumption Ppcypc also contribute to the HV 
battery power demand Pgattdem: 


PBatt.dem = Prrac.elec + Petc + Pac + Pbcipc (4.17) 


In case of recuperation Pratt dem can become negative indicating that power is 
transferred to the battery. 


The HV battery consists of four cell strings arranged in parallel. Each cell 
string is composed of 100 lithium nickel-manganese-cobalt oxide (LINMC) 
cells connected in series [207]. The cells are liquid cooled and can be heated 
using PTC elements. The battery can provide an electrical power of 300 kW at 
a nominal voltage of 370 V [14]. Its nominal capacity of 38.3 kWh enables a 
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vehicle range of more than 200 km [16]. The grid charging power Peria is up 
to 22 kW with AC via the battery charger inside the vehicle and up to 100 kW 
with DC [14]. 


According to [71] the terminal voltage Ury of a LiNMC battery can be de- 
scribed using the following model stated in [134]: 


Urv = Uocv(SOC) -R-I 


a [> 120 (4.18) 
Rt 1<0 


Hereby Uocv is the open-circuit voltage as a function of state of charge (SOC). 
The internal ohmic resistance R can be separated into charge resistance R* and 
discharge resistance R~. The battery current / is assumed positive for discharge 
and negative for charge. Uocy, R* and R™ are determined by cell testing [134] 
and have been provided for the HV battery of the research vehicle in look-up 
tables. 


Figure 4.4 shows the equivalent circuit from [134] and the approximate open- 
circuit voltage derived from [71]. A comprehensive literature review of battery 
models and their application to LINMC batteries is given in [71]. 


According to (4.18) the internal ohmic resistance of the battery causes a HV 
battery power loss PBatt,Loss, during a HV battery power demand Ppatt,dem: 


Shy) 
Uocv :I- R -I° , Ppat,dem > 0 
=:Pp, =: Pgatt.Loss >0 
Pgatt.dem = Ury - I = 4 "Prmu>0 Ga (4.19) 
Uocv -I- RI > PBatt.dem <0 


=: Par <0 =:PBatt,Loss >0 


This means that if the HV battery power Pgat is positive, i.e. the battery is 
discharged, the battery must provide PBatt,Loss in addition to Pgattdem in order 
to compensate for its ohmic losses. As a result, the SOC will decrease faster. 
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Figure 4.4: High voltage battery characteristics. Top: Equivalent circuit according to the "sim- 
ple model" in [134] with state of charge (SOC), open-circuit voltage Uocv, charge 
resistance R*, discharge resistance R7, battery current J and terminal voltage Ury. 
The diodes are assumed ideal. Bottom: Exemplary Ury of the LiNMC battery in the 
research vehicle depending on SOC and voltage across the internal ohmic resistance. 
Uocv was derived from [105]. 


Conversely, if Pgatt is negative because of recuperation, the SOC will increase 
slower because the absolute value of the charging power is reduced by Ppatt.Loss- 


4.4 Energy consumption and optimization 
approach 


O Coming from a temporal representation, the HV battery energy Epa needed 
for a route results from integrating Ppa Over time t, whereby tTrip denotes the 
trip time: 


rip 
E f Poan (t)dt (4.20) 
t=0 
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Adopted from [79]. One goal of the ALC is to compute an energy-efficient ve- 
locity trajectory. For reducing the energy consumption of a BEV three different 
cost functions in time domain are compared in [81]: 


First, motivated by reports in literature that a velocity profile with low accel- 
eration and deceleration is beneficial for energy savings, a penalization of the 
absolute value of acceleration is investigated. Due to the resulting reluctance 
to change velocity limits, this approach is suitable for the considered ACC ex- 
ample only to a limited extent. The approach offers the least energy savings in 
the comparison (8.4 % and 4.1 % depending on the scenario). This approach 
is also used in [107] and has the benefit that no power consumption model is 
needed. 


Second, the difference of Pray from its minimum value is penalized, which 
means that the recuperation maximum is seen as the optimal state. This ap- 
proach leads to unnecessary strong and frequent recuperations, that are fol- 
lowed by corresponding accelerations. As losses occur in both recuperation 
and traction mode, this approach is the second least effective (10.6 % and 5.2 % 
on average). 


The highest effectiveness (11.9 % and 5.6 % on average) comes from avoiding 
deviations of Pat from zero only for the case of traction and to not consider 
recuperation. 


O The rationale of the energy consumption optimization approach in this work 


is to avoid power losses between battery and wheels denoted by PL oss: 


PLoss = PBatt — Pirac,mech (4.21) 


The corresponding efficiency measure Piracmech/ Pau is known as energy effi- 
ciency of driving and tank to wheel efficiency [10]. 


PLoss can be expressed with respect to Pirac elec as in Figure 4.5. The data 
points result from computing the driving resistances and power train quantitites 
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Figure 4.5: Power loss between HV battery and wheels PLoss versus electrical traction power 
Pirac,elec for various combinations of vehicle velocity vyncı and vehicle longitudinal 
acceleration Vypelx- Adopted from [79]. 


for various combinations of vehicle velocity vyncı and vehicle longitudinal 
acceleration Vynel,x With road curvature and slope both equal zero. 


Poss increases with the absolute value of vyncı and Pyhelx. Arrows indicate how 
the operating points are shifted when vyncı and Pyncl,x are increased. The PLoss 
minima for each Pirac elec can be approximated by a parabola through the origin. 
The slope of the parabola increases with the absolute value of Pirac elec indicating 
that the efficiency of the power transmission decreases with increasing amount 
of transferred power. Adopted from [79]. 


Sensitivity analyses by [10] reveal that Prac mech/Ppatt is the parameter with the 
highest impact on the energy consumption in a BEV. Frequently, consumption 
models that require detailed knowledge of these usually only roughly known 
parameter values are used at the risk of inaccurate results [10]. To overcome 
this problem, this work applies data-driven models that learn such relevant 
parameters in aggregated form during vehicle operation. 
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O Due to lack of sensors PLoss is not known during vehicle operation and 
therefore cannot be used for a data-driven model. However, Purac elec can be 
derived from voltage and current sensor data of the power electronics and can 
serve as a proxy for Poss according to the above considerations. Therefore the 
energy optimization criterion in the cost function of the ALC is designed to 
penalize the absolute value of Pirac,eiec- Adopted from [79]. 


In comparison, the approach of this work comes closest to the third cited method 
but goes beyond in that it additionally penalizes the also imperfect recuperation 
mode. Furthermore, using Pirac.elec as Optimization criterion is assumed to be 
more effective than using acceleration and deceleration as the first cited method 
does. This is because for the same requested acceleration or deceleration the 
required Pirac.elec can differ significantly depending on road slope and vehicle 
velocity. 


O The power required for auxilaries such as light, ventilation, heating, cooling 


and radio can contribute to the total energy consumption for a trip significantly, 
in case of low velocity vyncı < 30 km/h. For vyne > 80 km/h the auxiliary 
power demand is negligible. It also strongly depends on environmental condi- 
tions and individual comfort preferences. Its 5 % and 95 % quantils are given 
by 0.2 kW and 1.3 kW [10]. Adopted from [79]. 


Furthermore, some components such as the AC compressor are only temporarily 
active and their activation can be scheduled. Therefore a predictive operating 
strategy that consideres both traction power demand and cooling power demand 
has much more energy-saving potential compared to for instance penalizing the 
sum of instantaneous AC power consumption and electrical traction power. A 
longitudinal control that includes operating strategies of auxiliaries is presented 
in [183]. 


O For these reasons this work focuses on optimization of traction consumption 


and neglects auxiliaries by assuming Pßatt,dem = Pirac,elec- Even without explicit 
consideration of power demand of auxiliaries, the penalization of Ptrac elec has a 
beneficial effect on efficiency because any additional consumption of auxiliaries 
further increases the HV battery power loss. 
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The goal is to penalize Ptrac elec only to such an extent, that the ALC avoids 
inefficient power peaks that have negligible effect on the trip time. Adopted 
from [79]. With large penalty for Pirac elec, the vehicle will drive very slow or 
even come to a standstill because of lack of Pirac,mech that is needed to overcome 
the driving resistance. Then the resulting trip time will strongly increase as 
well as the energy consumption because of the power demand of auxiliaries. 
However such large penalty for Ptrac elec is far beyond the driver acceptance. 


4.5 Vehicle models for simulation 
environment 


The preceding sections investigated the driving resistances that result from a cer- 
tain vehicle velocity vyhe1, vehicle longitudinal acceleration Vyhe1 x, road slope 
angle œ and road curvature « and stated how the corresponding HV battery 
power Par can be derived from the mechanical traction power Pirac,mech- The 
upper diagram in Figure 4.6 depicts this approach, which is referred to as back- 
ward simulation. Starting from the desired driving state backward simulation 
computes the state of components such as gear box against the direction of the 
effect chain. 


Backward simulation assumes that the power train can provide the power 
needed for the given driving situation. Due to the missing consideration of 
physical root causes, this approach is also denoted non-causal modelling [45]. 
Non-causal models include no control loop, therefore no controller is needed. 
However, backward simulations are only suitable for quasi-stationary simula- 
tions that do not take into account transient effects. Quasi-stationary states 
of components are usually derived by interpolation between values of look-up 
tables as mentioned in Section 4.3. 


Backward simulation will be applied for computing the resulting energy con- 
sumption for planned trajectories in Chapter 5 and for recorded drives in Chap- 
ter 6. The model will be referred to as open-loop reference model. 
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Backward simulation 


&, K, VVhch VVhcl,x Driving P trac,mech 
_— 
resistances 


Forward simulation 


Driving P trac,mech 


resistances 


VVhcl» VVhel,x 


Vdes» Ades Controller 


Figure 4.6: Sequences of forward (top) and backward simulation (bottom). 


The equations from the previous sections can also be rearranged for a forward 
simulation. The lower diagram in Figure 4.6 depicts this case. Forward simula- 
tion computes component states starting from the cause in accordance with the 
cause-and-effect chain. Therefore this approach is denoted causal modelling. 


Forward simulations enable dynamic simulations that incorporate dynamic com- 
ponent behavior, e.g. transient effects in the motors and battery. Components 
can be described using differential equations which are connected to each other 
according to the cause-and-effect chain. In this work the cause is the desired 
motor torque Tyes While vyncı,x is at the end of the cause-and-effect chain. vyncı,x 
is derived by rearranging (4.9) and used for computing vvncı for the next sim- 
ulation time step. Causal models include a closed control loop and therefore 
require a controller. 


Forward simulations will be used for testing the interaction between ALC and 
vehicle. The model will be referred to as closed-loop reference model. The 
controller inside the ALC chooses Taes such that deviations of vyncı and Vyncix 
from the desired velocity vdes and the desired acceleration Ages, respectively, 
ideally vanish. Vaes and Ages are specified by the planned trajectory of the ALC 
[48], [62, pp. 37-41]. 
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Figure 4.7: Course of the Weissach route (WR) that is used as reference route for the evaluation 


of the automated longitudinal control. The dotted arrow indicates the driving direction 
and the straight line the start and end point of the WR. 


4.6 Model of reference route 


O The reference route is a roughly 23 km long circuit around the village Weis- 
sach in Southwest Germany, hence denoted Weissach route (WR). Figure 4.7 
depicts on a map the WR, that comprises urban sections as well as country 
roads. The legal speed limit varies between 30 km/h and 100 km/h and the road 
slope ranges from -8 % to 10 %. Adopted from [79]. 


Figure 4.8 displays the corresponding map data consisting of legal speed limit, 
road curvature and road slope as well as the therewith derived elevation profile. 
The map data was extracted from a navigation system and will be used in back- 
ward and forward simulations with the closed-loop and open-loop reference 
model, respectively. 
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Figure 4.8: Route data of Weissach route (WR) consisting of legal speed limit, road curvature and 


road slope used for simulation as well as the therewith derived elevation profile. The 
position is measured from the start of the WR (c.f. Figure 4.7). 
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4.7 Vehicle models for automated 
longitudinal control 


O This section presents two models for use within the ALC system. For a given 
driving situation, the first model provides the required traction force Frac and 
the second one the required electrical traction power Pirac,elec- Both quantities 
depend on the driving resistances, which in turn are affected by the given driving 
situation as well as parameters. Adopted from [79]. 


According to sensitivity analyses in [10], the second most relevant parameter 
for the accuracy of a vehicle model is fr. In the comparison, the model is only 
affected by p for vyncı > 100 km/h, whereas cw and A are negligible in the 
investigated situations. myncı is only relevant for hilly trips and trips that lead 
to higher elevation. 


O During vehicle operation these parameters are not exactly known. Myhel 
changes with additional passengers or luggage, cy and A when the convertible 
top is being opened and f, with tire temperature. 


Therefore adaptive models are applied, which estimate these parameters during 
vehicle operation either explicitly or in aggregated form. In order to be able 
to adapt a model, its inputs and outputs must be quantities that can be derived 
from signals on the Controller Area Network (CAN) bus. Adopted from 
[79]. However, only few quantities of the power train that are relevant for 
this purpose, are measured. Therefore the adaptive models represent the power 
train properties only on an aggregated level and cannot state losses in individual 
components, e.g. gear box losses. 


Compared to the vehicle reference models for the simulation environment both 
vehicle models for ALC are simplified in that they first do not take into account 
road curvature explicitly and in that they second do not distinguish between 
vehicle mass and rotational masses of the power train as (4.9) does. 
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Some publications take into account the inertia of rotating parts of the power 
train [96, 188] for more accuracy but the majority does not [12, 91, 102, 127, 
195, 204]. Due to lack of a fly wheel, less complicated gear boxes and smaller 
motors, the inertia of rotating parts are usually less in a BEV than in a conven- 
tionally driven vehicle. As the inertia of rotating parts is usually not known and 
elaborate to determine, publications that consider it, mostly increase the vehicle 
mass by 5%. In the research vehicle, the equivalent masses add up to roughly 
5.6 % of the vehicle mass. Varying the mass factor between 0 % and 5 % had 
almost no impact on the results in the sensitivity analyses for a BEV in [10]. 


O The simplification regarding rotational masses allows using the sensor lon- 
gitudinal acceleration a, as a model input. a, is measured by the acceleration 
sensor and is influenced by both change of velocity and road slope angle: 


ax = Vvhelx + 8: sin(a) (4.22) 


Adopted from [79]. The required traction force and the required electrical 
traction power can differ for the same a, even if all other influence factors are 
kept constant. This is because in the case that a, results from driving at a slope 
with constant velocity, no rotational acceleration of the power train occurs but 
it does if ax results from change of velocity at zero slope. The adaptive models 
cannot distinguish these cases. The simplification using ax is frequently used 
for driving resistances parameter estimation [140, 170, 178]. 


A comprehensive review of vehicle energy consumption models including an 
analysis of influence factors on the energy consumption is given in [118]. Fur- 
thermore, [118] classifies consumption models into white-box, gray-box and 
black-box models. White-box models are based on high knowledge of the 
underlying system and incorporate detailed descriptions of subsystems. In con- 
trast, black-box models only learn an input-output pattern from provided data 
and require no system knowlege. Gray-box models partly use system knowlege 
and are partly driven by data. 
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Based on the inputs to a black-box consumption model, [118] distinguishes 
engine-based consumption models with engine torque and engine speed as 
inputs, vehicle-based consumption models with instantaneous vehicle speed 
and acceleration as inputs and modal-based consumption models which use 
operating modes such as idling, accelerating or cruising as inputs. 


Based on the criteria in [118], the proposed traction force model can be regarded 
as a gray-box model and the proposed model for the electrical traction power 
as a vehicle-based black-box consumption model. 


Similar models can be found in literature. For example, [48] models the tractive 
force without cornering resistance and the mechanical traction power of aBEV 
similar to (4.11) and (4.12) depending on instantaneous speed, acceleration and 
road grade information. The driving resistance parameters are assumed to be 
known and constant. A parameter that describes the efficiency of recuperation 
is estimated with the LS method. 


In contrast, [51] takes a more aggregated approach and models the traction force 
with a quadratic polynomial as a function of vehicle velocity. This traction 
force model is then used for a multivariate model of power demand and energy 
consumption of a BEV that depends on vehicle velocity and acceleration. In 
[194] the polynomial regression model additionally uses the SOC as an input. 
Two KAF algorithms are compared to RLS for representing Pirac,elec depending 
on a, and velocity in [146]. 


4.7.1 Adaptive traction force model 


O The ATFM answers the question how much traction force is needed in order 
to fulfill the driving demand of the longitudinal control in the current situation. 
The ATFM is used as a pilot control in the controller depicted in the lower 
diagram in Figure 4.6. The controller computes a motor torque demand such 
that the vehicle tracks the planned velocity trajectory. The ATFM is based on 
the following simplification of the longitudinal traction force equation (4.11) 
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that neglects the cornering resistance (4.6) and the equivalent masses of the 
power train in the inertial force (4.9): 


Frac = Fion + Finert,simplified + Fa + Fair 


Mvyhcl * 8 * fr * COS(@) + Myhel * Vvhel,x (4.23) 


: P 2 
+ Mynci + 8 > Sin(@) + 7 - A: Vvhel 


For small slope angles cos(@) ~ 1 applies. Therefore the first summand can 
be approximated by the rolling resistance constant Fo when additionally the 
vehicle mass in this summand is assumed to change only slowly. Moreover, the 
second and third summand are merged using (4.22). The traction force Frac iS 
not measured but it can approximately be computed as 


ig 


+ TEM.RA . (4.24) 


Firac = TEM,FA * , 
Fdyn,FA Tdyn,RA 


whereby the braking torques in the power train are neglected in contrast to 
(4.16). Tema and Temra are no measured quantities either. These signals are 
computed by the motor ECU from the measured motor voltages and currents 
using look up tables. 


With these adaptions (4.23) in matrix form reads 


Frac = (1, ax, Vyhel”) - (Fo, Mynci» EA)", (4.25) 
ee 2 


U 


=:C Vhel BE 
i IC: 


whereby Cynci is the vehicle motion vector. The vehicle parameter vector Xx yhel 
summarizes the driving resistance parameters. These parameters are generally 
not fully known. Therefore x ypc, needs to be estimated such that the estimation 
minimizes the residual between the traction force computed with (4.24) and the 
traction force according to the model output from (4.25). Adopted from [79]. 
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This optimization problem is solved with a KF (c.f. Subsection 3.3.2) and the 
following model equations: 


Xvhel,p = XVhel,p-1 + @p (State equation) (4.26) 


Firac,p = Cvhel * XVhel,p + Up (Measurement equation) (4.27) 


p denotes the time step, w the process noise and v the measurement noise. 
Frequently, algorithms such as RLS, KF, EKF are used to determine the driving 
resistance parameters [140, 170, 178]. O For the ATFM a Stenlund-Gustafsson 
M-Kalman filter described in [145] is applied. It uses a regularization from 
Stenlund and Gustafsson [168]. This regularization method sets lower bounds 
for the covariance matrix in such a way that the estimated parameters are kept 
constant in phases of low excitation, i.e. when Firac and C ype temporally remain 
constant, such as while driving with roughly constant velocity and road slope. 
Low excitation also applies to vehicle standstill but for this phase a simple 
criterion for pausing the adaption can be defined. Adopted from [79]. The 
ATFM was created by F. Bleimund [22, pp. 24-29] within e-generation and 
reused by the author of this work without modifications. 


Literature also proposes methods that focus on fusing the information of differ- 
ent sensors [74, 198] as well as methods that divide the estimation into several 
stages [101, 199]. 


4.7.2 Adaptive electrical power model 
Model interfaces and features 


O The AEPM answers the question how much electrical traction power Pirac,elec 
is required to fulfill the driving demand. The ALC uses the AEPM during 
trajectory optimization in order to derive trajectories that require few electrical 
traction power. As shown in Section 4.4, this corresponds to a low power loss 
between HV battery and wheels. 
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Prracelec can be determined from the measured voltages and currents of the 
electric motors that are available on the CAN bus. Apart from the power losses 
in the drive train, Ptrac elec equals the mechanical traction power Pirac,mech Which 
can be computed as the product of traction force Firac and vehicle velocity vvncı 
(c.f. (4.12)). According to Subsection 4.7.1, Fira can be modelled as a function 
of the CAN signals vyncı and ax. Therefore the AEPM relates vypcı and ax to 


P trac,elec: 


Pıepm = AEPM (vyncı, ax) (4.28) 


Features of power consumption representation 


A substantial fraction of Pirac.mech 18 needed to overcome the acceleration resis- 
tance Myhcl * VVhcl,x from (4.9). The corresponding acceleration power myncı * 
VVhcl,x*VVhcl depends on an unseparable product of Vvhe] x and vvncı. For negligi- 
ble road slope angle @ these quantities are similar to the model inputs because of 
(4.22) (c.f. [51]). Therefore the linear model structure used in Subsection 4.7.1 
for the ATFM is not suitable for the AEPM. Adopted from [79]. 


The power consumption representation in time domain is highly nonconvex, 
as stated at the example of a Pra = f(Firac, Vvhel) representation in [81]. An 
approximation of Ppat with a convex function can be used as a cost function in 
an convex optimization problem that can be solved more easily than a nonlinear 
optimization problem. The quadratic approximation accuracy for various repre- 
sentations is also investigated in [81]. Due to the less prominent nonconvexity 
of the representation in space domain, the highest accuracy is observed for the 
(PBatt/VVncı) = f (Frac, Vyhl”) representation, which is the energy consumption 
in Watt seconds per meter. In space domain, however, vehicle dynamics and 
the vehicle control problem are more nonlinear [81]. 
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Algorithm for model adaption 


O The AEPM is based on the Fixed-Budget KRLS (FB-KRLS) algorithm de- 
scribed in [180], which is a representative of KAFs introduced in Section 2.4. 
The source code of FB-KRLS is available at [181]. With each new data point, 
FB-KRLS adds a new kernel support vector in the vvpcl — ax plane, updates 
the approximation and discards the least important kernel support vector to 
keep the required memory constant. During model evaluation the model output 
Parpm results from computing the sum of all kernels at (vvncl, dx), weighted 
with their corresponding control points. Adopted from [79]. 


The model captures aggregated power train losses that are constant or depend 
ON Vvnel Or ax. If vehicle parameters or power train properties change, the 
AEPM will adapt itself accordingly. 


While [81] approximates the nonlinear structure by a convex function to get 
a convex optimization problem, a KAF approximates the nonlinear structure 
using a nonlinear function and still offers a convex optimization problem for 
model adaption. For this kind of convex problems there are solvers that con- 
verge in a defined time period, which is important for real-time model adaption. 


RLS and two fixed-budget KAFs, called KRLS Tracker (KRLS-T) and QKLMS- 
FB, are applied to (4.28) in a very similar approach in [146]. KRLS-T achieves 
the highest accuracy, followed by QKLMS-FB. As expected, the linear RLS, 
proves not suitable for strong nonlinearity and achieves the lowest accuracy. 


All three mentioned KAF use the Gaussian kernel defined in (2.21). The KRLS- 
T algorithm published in [179] in 2012 is a KAF algorithm for time-varying 
regression that includes a forgetting factor that can handle non-stationary scenar- 
ios. In contrast, the FB-KRLS published in 2010 stems from a sliding-window 
approach with the improvement that it maintains not the kernels for the latest 
but for the most important data. Its forgetting mechanism is not optimized 
towards non-stationary scenarios [181]. 
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The idea to model the power train characteristics with a FB-KRLS originates 
from S. Rhode [146], who also provided a first script in MATLAB that performs 
this task. Based on this, both researchers continued the investigations inde- 
pendently from each other, so details differ. The content of the remainder of 
this subsection, except from the normalization described below, stems from the 
author of this work. The novelity of this work regarding FB-KRLS is limited 
to the application of this kind of black-box model for trajectory optimization 
and within the ALC along with required problem-specific adaptions such as 
the data shift as well as the determination of model hyperparameters, both also 
mentioned below. 


Data transformation 


O For various (Vvncl, 4x) combinations the trajectory optimization algorithm 
evaluates the AEPM and takes decisions based on the AEPM outputs. The 
standard FB-KRLS output is zero in border areas where no data points have 
occured yet. This means that if the trajectory optimization evaluates the model 
for high velocities or accelerations that are beyond the capabilities of the vehicle, 
the AEPM states Pracelec ~ 0. To avoid that results of evaluations of the model 
in its border areas appear as efficient and therefore favorable driving states to 
the trajectory algorithm, the model outputs must have large absolute values 
in border areas. This can be achieved by shifting Pracelec in training data by 
subtracting a large offset value. 


After the shift, all components of training data (Vynel» 4x, Ptracelec) and evalu- 
ation data (Vvnel, ax) are additionally scaled to a ranges between -1 and 1 by 
multiplying each of them with a constant between O and 1. This process is 
called normalization. Frequently better aproximation quality is observed for 
models that work with normalized data. 


Shift and normalization form the transformation. Results of model evaluations 
are retransformed by a denormalization and backshift operation before they are 
provided to the trajectory optimization algorithm. 
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Figure 4.9: Electrical power Pagpm according to the adaptive electrical power model (AEPM) 
depending on sensor longitudinal acceleration ax for various vehicle velocities vyhe]. 
The gray shaded area is the ALC operating area, max Pirac elec the maximum traction 
power and min Pirac,elec the maximum recuperation power. Adopted from [79]. 


Figure 4.9 depicts the retransformed model output as a function of a, for various 
nonnegative vehicle velocities. Negative velocities are not shown for simplicity 
and because with ALC the vehicle only drives forward. For medium negative 
ax the Pirac elec 18 negative indicating that power is recuperated into the battery. 
The gray shaded area indicates possible operating points of the vehicle while 
driving with activated ALC under the assumption of no slope. As the ALC 
limits are defined with respect to Vvnelx, the width of the area will change in 
presence of slope according to (4.22). 


In border areas that exceed the vehicle capabilities the shift leads to very large 
AEPM outputs of up to three times the maximum Pirac elec. AS a result, the tra- 
jectory optimization algorithm avoids trajectories that include such unreachable 
(VVhel, dx) combinations. Adopted from [79]. 


In case of negative slope peaks or strong braking, the simple shift leads to a 
zero-crossing Of Pirac.elec at roughly —7 m/ s? < ax < —5 m/s? before Pirac elec 
increases to large values. This case can still mislead the trajectory optimization 
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algorithm but for the given setup it did not present a problem. As Figure 4.8 
indicates, on the WR peaks in the course of the road slope are stronger in the 
positive direction than in the negative direction. This asymmetrical road slope 
distribution causes that in AEPM evaluations the maximum absolute value of 
a, is greater for a, > O than it is for a, < O and that border areas on the left 
hand side of Figure 4.9 are not reached. 


To prevent that negative slope peaks cause issues on other hilly routes, the 
shift should be improved by subtracting a plane from the model input data 
instead of a constant. This plane can be designed so that it goes through the 
origin and is tilted upwards in the directions of the vyhe and ax axes, e.g. 
—C1 * VVhel — C2 * Ax + €3 * Pirac,elee = O with constants c1, c2, c3 > 0. By choosing 
these constants appropriately, physically unreasonable zero-crossings of PAEPM 
can be avoided. 


Determination of model hyperparameters 


The question is how kernel variance 07? in (2.21) and kernel count M should 
be chosen for the AEPM. For better differentiation from other parameters like 
control points, da? and M are called hyperparameters [146] and determining 
them is part of the model selection process [72, pp. 128-130]. A model with 
many parameters usually provides more accurate results at the cost of increased 
effort for parameterization and evaluation. 


Model selection criteria, also called information criteria, assess models and 
thereby can help to find a trade-off between model complexity and model qual- 
ity. The general structure of an information criterion is: 


SSE 
criterion = P + In (55) + penalty term (4.29) 
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In(-) is the natural logarithm, P the number of observations and SSE is the sum 
of squared errors with 


$ 
SSE = Y (yp- 9p) - (4.30) 
p=1 


y denotes the measurement and $ the model output. The selected criterion is 
computed for a set of trained models and the model with lowest criterion value 
is chosen. Several information criteria have been introduced, which are based 
on different principles and differ in the penalty term. In a comparison in [89] 
the Hannan-Quinn information criterion (HQC) performs best: 


HQC =P - In (=) + M - In (In(P)) (4.31) 


Number of model parameters M and kernel variance o* from (2.21) span the 
hyperparameter space for the search for the most suitable AEPM. For each 
combination of M = 10,20,..., 100 and 0? = 1/1, 1/2, . . . , 1/30 a model is trained 
and tested to determine the model with the lowest HQC. 


Data for model training and test was created using the open-loop reference 
model. Thereby the vehicle longitudinal acceleration was sampled with vyncıx = 
-10,-9.8,..., 10 m/s” and the vehicle velocity with vyna] = 0,0.5,...,56 m/s. 
Slope and curvature were assumed equal zero. Operating points that are be- 
yond the vehicle capabilities according to the open-loop reference model were 
removed from the data set. The data set that comprises all possible vehicle 
operating points (PVOP) in random order is split into a PVOP training data set 
and a PVOP test data set. The randomly chosen PVOP test data set comprises 
10 % of the original PVOP data set and the remaining data points form the 
PVOP training data set. 


Additionally, a WR data set was created from the closed-loop reference model 
output during a simulated drive with the ALC on the WR including slope and 
road curvature information. The WR data set was split analogously into WR 
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training data set and WR test data set. PVOP data set and WR data set contain 
roughly the same amount of data points. 


The upper diagram of Figure 4.10 depicts the best HQC as well as the corre- 
sponding o for each investigated kernel count. Only for models that are both 
trained and tested using data from the PVOP data set the HQC course shows 
a local minimum at 50 kernels along with comparatively large kernel variance 
values. 


The lower diagram of Figure 4.10 depicts the corresponding 95 % quantile Q95% 
of absolute values of the relative error ere; between the output of the AEPM and 
the reference model for PVOP and WR, respectively. With increasing M the 
95 % quantiles decrease. 


Models that are both trained and tested with PVOP data achieve better HQC 
values and lower quantiles than models that are trained with PVOP data and 
tested with WR data. The PVOP data set contains much more dynamic driv- 
ing states than the WR. Furthermore, only the PVOP test data set evaluates a 
model that has been trained with the PVOP training data set in its border areas. 
Therefore, the worse result for the WR test data set seems counterintuitive at 
first. 


However, only the WR data set considers road curvature and road slope. On 
the WR the road slope is up to 10 %. The AEPM does not use road curvature 
input and summarizes road grade and change of vehicle velocity in the input ax. 
Hence, the results indicate that road geometry has more influence on the quality 
of the model output than the driving style. However, for M > 40 the quantils 
are below 2 % for all considered combinations of training and test data sets. 
Therefore not considering road curvature in the AEPM explicitly is acceptable. 


Since the goal is not to optimize the model parametrization with respect to 
the WR but to find a suitable model for the vehicle in general, the best model 
parametrization according to HQC for the PVOP test set was chosen. This is 
a model with M = 50 and 0? = 1/9 that achieves O95%(lereıl) = 0.94 Yo for 
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Figure 4.10: Results of full factorial search for best model hyperparameters kernel count M and 
kernel variance 0”. Top: Both the best HQC and the corresponding o? are plotted 
versus M for each combination of training and test using the data sets from possible 
vehicle operating points (PVOP) and Weissach route (WR) . Bottom: 95 % quantile 
Q95% of absolute values of relative error ere; between output of AEPM and reference 
model versus M. 


PVOP and Oosx(lereıl) = 1.5 % for WR. This AEPM model is used in the 
remainder of this work and depicted in Figure 4.9. 


With the hyperparameter optimization technique from [146] the KRLS-T accu- 
racy increases for M < 50 and converges for more kernels. Therefore M = 50 
is considered the optimal KRLS-T hyperparameter, analogously to this work. 
However, with 0? = 6.25 the reported kernel variance is different. Reasons can 
be contrasting data sets and differences in the KAF algorithms themselves. 
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This chapter presents a trajectory planning approach based on the B-spline ap- 
proximation methods from Chapter 3. Section 5.1 states how an upper speed 
limit for the route ahead is created from map data. Section 5.2 includes con- 
siderations on how to represent the planned vehicle velocity and Section 5.3 
describes the trajectory optimization. Section 5.4 proposes an extension of the 
trajectory optimization problem for taking into account the required electrical 
traction power. Section 5.5 states ways to enforce constraints on the trajectory 
beginning and Section 5.6 summarizes the scientific contribution of this work 
regarding trajectory planning. 


5.1 Generation of upper speed limit 


An upper speed limit for each meter of the road section ahead of the vehicle 
is computed using map data. The upper speed limit takes into account legal 
restrictions, limitations from driving dynamics as well as safety and comfort 
requirements and serves as input data for the trajectory optimization process. 


The map data consists of three P x 1 vectors, each with component index 
p = 1,...,P. The vectors state the meter-discrete courses of road slope y, 
road curvature x and legal speed limit VLim Law for the road section ahead of the 
vehicle in driving direction. The first component of each vector refers to the 
current vehicle position. P is called the length of the electronic horizon. 


Figure 5.1 depicts various speed limits that occur during generation of the upper 
speed limit and are described in the following paragraphs. 
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VLim,Law VLim,Speedo VLim,Curve 
VLim,Map VLim,Map,v 


VLim,Crest 
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Figure 5.1: Computation of upper speed limit with respect to position measured from vehicle at 
km 11 of Weissach route in driving direction. Corresponding route data depicted in 
Figure 4.8 between km 11 and km 13. Legal speed limit vLim,Law, legal speed limit in- 
cluding speedometer offset vLim,Speedo» Speed limit resulting from curvature Vy im,Curve 
and speed limit resulting from crests VLim Crest are intermediate quantities. The speed 
limit resulting from map data VL im Map is the minimum of VLjm,Speedo» VLim.Curve and 


VLim.Crest- End result is the speed limit from map data with desired acceleration 
VLim,Map,v - 


The course of the legal speed limit vLim,Law forms the basis of the upper speed 
limit. The speedometer in a vehicle includes an offset and therefore indicates 
a velocity that is slightly higher than the actual vehicle velocity. For better 
acceptance by the driver, the legal speed limit including speedometer offset 
VLim,Speedo 18 used instead Of vLim,Law. This causes that the speedometer will 
indicate roughly the same velocity value as the traffic sign when there are no 
relevant restrictions other than VLim Law present. 


Due to tight curves, crests or comfort requirements the upper speed limit needs 
to be corrected starting from vLim,Speedo further downwards. 


116 


5.1 Generation of upper speed limit 


Driving with vehicle velocity vyhe on a road with curvature «x causes a lateral 


Vo - K (c.f. (4.5)). A maximum absolute value of lateral 


acceleration Vvyhel y max 18 specified, which leads to the speed limit resulting from 


acceleration Vynely = 


curvature VLim,Curve given by 


VVhel,y,max 
VLim,Curve S —— e (5.1) 
[kl 
Thereby | - | denotes the absolute value. A comfortable lateral acceleration 


maximum varies with vehicle speed [153]. As a first draft a lateral acceleration 
table Vynel y max (V) was created using a characteristic curve from [153]. In test 
drives the characteristic curve was further adapted to suit the vehicle charac- 
teristics. The lateral acceleration table consists of (Vvncl,y.max» v) supporting 
points because these quantities can be specified more conveniently than x in 
test drives. With (4.5) the corresponding x values are computed. Thereafter 
Ývhcl,y,max for a given « is determined from map data by linear interpolation and 
VLim,Curve 18 derived using (5.1). 


Crests on the route cause another restriction on the upper speed limit. From road 
slope data an elevation profile is obtained and it is calculated how far ahead the 
road can be seen according to geometrical considerations. It needs to be ensured 
that the vehicle can always come to a standstill within the sighting distance 
ASsighting. Assuming that the absolute value of the maximum deceleration is the 
gravitational constant g leads to the speed limit resulting from crests Vim Crest 


given by VLim.Crest < 42 - ÁSSighting * 8- 


The speed limit resulting from map data vLim.Map is the minimum of VLim,Speedo» 
VLim,Curve AMA VL im. Crest- Usually VLim,Crest 18 high compared to the other quanti- 
ties and therefore rarely determines VLim,Map- 
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In the last step the desired longitudinal acceleration ages is taken into account. 
The desired acceleration and deceleration values are derived by interpolat- 
ing with respect to velocity v in separate acceleration tables for acceleration 
Ades pos(V) and deceleration ddes,neg(v) < O. Enforcing the constraints 


Vp+1 S a/Vp +2-As- Adespos(Vp), P = 1,2,...,P—1 


vp-1 < „vp 2: As - daesneg(Yp), p = P,P- 1,...,2 


(5.2) 


with As = 1 on VLimMap gives the speed limit from map data with desired 
acceleration VLim,Map,v, Which is the output of the whole procedure. 


The content described above in this section originates from F. Bleimund [22, 
pp. 29, 30], who implemented it during e-generation. The adaptions done 
during e-volution are limited to enhancements of the acceleration tables and 
adaptions of their entries during test drives. 


Other works such as [114, 139] additionally specify a lower speed limit as a 
certain percentage of the upper speed limit. Both limits then form a so-called 
driving tube that defines the solution space of valid trajectories. An alternative 
term is driving envelope [44]. This work specifies no explicit lower limit but 
enforces a nonnegativity constraint on velocity trajectories as Section 5.5 will 
state. 


5.2 Representation of vehicle velocity in time 
domain 


There are several ways to define a course of velocity that stays below the upper 
speed limit. Figure 5.1 illustrates that in general only a functional representation 
with many degrees of freedom can sufficiently adapt to the upper speed limit 
over several kilometers. 
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Section 2.3 investigated the properties of several function types and the advan- 
tages of spline functions over polynomials. Based on the conclusions drawn 
there, the velocity trajectory will be represented using a cubic B-spline function. 
This function is twice continuously differentiable which means the trajectory 
velocity vryy is smooth. This is beneficial for driving comfort but not sufficient 
for jerk-free driving behavior. If the B-spline function has many degrees of 
freedom in a short time period and unfavorably chosen control points, it can 
lead to an oscillating and uncomfortable velocity course. 


A further helpful feature is that in each interval a B-spline function lies within 
the convex hull of the control points that are relevant for this interval (c.f. 
Section 2.3 and Section 3.1). This feature allows to enforce vrsy > O by simply 
restricting the control points to nonnegative values. 


Since the upper speed limit depends on the position, a straight-forward approach 
is to define the trajectory with respect to the position as well. This was done in 
[22, 139]. 


RBA described in Subsection 3.3.3 allows to adapt the control points of a B- 
spline function iteratively such that it approximates a set of data points in the 
WLS sense. With RBA a spatial velocity trajectory vryy(s) that depends on the 
position s can be adapted to the upper speed limit by solving the LS problem 


P 


A ; 2 
x= argmin X, (VLimtıy.p = vrıy (Sp)) (5.3) 
x = 
using the data set 
(Sp VLimTIY,p) Sp =P 1, p=1,2,...,P, (5.4) 


in which Vy im,Tyy,p 18 derived from vLim,Map,y as (5.12) will define. A spatially 
defined trajectory can easily be compared against the upper speed limit and 
adjusted downwards at any position without risking to violate the upper speed 
limit at another position. However, the spatial trajectory definition has at least 
four flaws: 
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First, at low speeds the position s changes slowly in a given time interval. If 
the knots of vryy(s) are spatially equidistant, the trajectory can only represent 
a comparatively reluctant driving behavior at low speeds, whereas at highway 
speeds it has an excessive number of degrees of freedom in each time interval 
that facilitates undesired velocity oscillations and increases the computational 
effort unnecessarily. 


RBA only chooses the control points and requires knots as an input. Regarding 
RBA this means that simply specifying spatially equidistant knots is not suitable. 
This was also confirmed by real test drives with a spatial trajectory planning 
approach based on RBA. Therefore a knot placement procedure was added that 
computed spatial knot positions that were equidistant with respect to time. The 
degrees of freedom were then spatially dense at low speeds and with increasing 
vryy their density decreased. 


There are also solutions for adapting both knots and control points when fitting 
a B-spline function to data [36, 93, 130] but this problem is a nonconvex one 
with many local minima and therefore difficult to solve [15]. For example, 
[158] solves a reduced nonlinear B-spline LS approximation problem, in which 
only the knots are optimization variables using a generalized Gauss-Newton 
method, whereas [141] applies the LM algorithm from Subsection 3.4.1. A 
review of state of the art methods for this kind of problem is provided by [43]. 


Second, a spatial velocity trajectory cannot represent driving off again after 
coming to a standstill because if vryy(s*) = 0, at any position s*, then s* cannot 
be left anymore. Hence, starting from a standstill requires a different trajectory 
planning approach. 


Third, the time ¢* at which a vrsy(s) trajectory reaches a certain position s* 
requires to compute the integral t(s*) = f: 2 ds. For a cubic function 
vrjy(s) the solution can be determined using a partial fraction decomposition 
(PFD). Since the solution of PFD depends on the control point values, a system 


of equations needs to be solved repeatedly. 
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Fourth, derivates of vryy(s) are no common quantities with known interpreta- 
tion. Common quantities refer to time and deriving them from vryy (s) is costly. 
For example, the derivative vryy’(s) is not the trajectory acceleration arıy (s) 
but the pseudo trajectory acceleration [55]. atyy needs to be computed from 
vryy (s) as a product: 


dvriy(s) 3 ds dvriy(s) 


Er Er ÓN vryy(s) : vryy (s) (5.5) 


arıy(s) = 


Hence, aryy cannot be repesented by the linear combination of basis functions 
and control points from (2.17) and (3.3), respectively, anymore. Equation (5.5) 
gives oscillations in the acceleration course, which become larger as the knot 
distance is decreased to allow for more agile vehicle behavior. The reason is that 
the result of (5.5) is no B-spline product and therefore also does not fulfill the 
convex hull property. Ways to correctly compute product functions of B-spline 
functions are given in [121, 133]. The approach in [133] relies on intermediate 
transformations to the Bézier functions (c.f. Section 2.3). A sliding window 
B-spline multiplication algorithm is presented in [30]. 


A time dependent velocity trajectory overcomes these flaws. First, the knots 
can simply be chosen temporally equidistant. Then the trajectory represents 
a vehicle behavior whose agility does only depend on the chosen constant 
knot density but not on vehicle velocity. Second, driving off from a standstill 
poses no problem because ¢ still changes during a standstill and so can vryy (t). 
Third, the corresponding trajectory position sryy (t) is obtained efficiently by 
integration: 


eee i uta: (5.6) 


Fourth, meaningful quantities like trajectory acceleration and trajectory jerk 
arise from a linear combination of basis functions and control points. 


For these reasons a time dependent trajectory vrjy(t) is chosen. The problem 
of adapting vrsy(f) to the upper speed limit is formulated analogously to (5.3) 
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as an approximation problem but an additional change of representation space 
from position to time is included: 


P 
x = arg min) (veim, (stiv(tp)) = vp) (5.7) 


x p=1 
The trajectory position stjy is given by (5.6) and computed as in (3.10). 


In contrast to (5.3), the error vLim,try (Stry) — vrry in (5.7) does not linearly 
depend on the control point vector x. Instead, x is linked to the error by the 
nonlinear upper speed limit, which makes (5.7) a NWLS problem. 


The described approach is similar to parametric B-spline curve approximation. 
As mentioned in Section 2.3, a B-spline curve is a generalization of a B-spline 
function that uses multi-dimensional control points. These enable an indepen- 
dent curve shape in each dimension such that the curve can represent edges of 
a geometrical object [3]. The velocity trajectory approximates (Sp, VLim,TIY,p) 
data points with respect to the parameter t. However, sqjy and vryy are coupled 
via (5.6) and by still using a B-spline function with scalar control points this 
constraint can be enforced. 


A disadvantage of the temporal trajectory definition is that road grade, speed 
limit and road curvature are position dependent, whereas the position of the 
vehicle depends on the future course of velocity. Therefore a velocity optimiza- 
tion requires a vehicle position estimate [82]. Specifically for the approach in 
(5.7) this means that if any component of the solution £ of the optimization 
problem is modified, the trajectory will likely exceed vi imap,» at some later 
point in time. In contrast, with vryy(s) control point values can always be 
reduced without risking to violate vLim.Map,» afterwards. 


For completeness a styy(t) trajectory is considered as well. Against such a tra- 
jectory argues the fact that stıy is unbounded and may take large values for long 
trajectories. This encourages numerical problems and instable computations in 
optimization algorithms. 
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5.3 Trajectory optimization 


o Trajectory optimization involves adapting the vryy(t) trajectory function to 


the data set in (5.4) that is created using the upper speed limit vLim,Map,v- VIIY (£) 
is defined according to (3.3) with degree d = 3, knot vector k and control point 


vector x. K given by 


K = (Kj, K2,...,KK) 


(5.8) 
= (At, : (-d), At. (=d+1),..., Atk: (=d+K-1)) 


has equidistant and strictly monotonously increasing entries. At, denotes the 
constant temporal distance of neighboring knots. Due to the choice of x, the 
trajectory can be evaluated for t > 0. 


t is discretized by a constant temporal distance of neighboring data points Aty: 
tp =(p-1)-Ath p=1,...,P (5.9) 


The approximation problem that is solved reads 


P 
Í = arg min ) (e; [vscı (stiy(tp)) = wre) 
poe (5.10) 


a 2 Y. 2 
+Rz! > (asa - anıy(tp)) +R! > (inv(tp)) }. 
arıy denotes the trajectory acceleration and jryy the trajectory jerk. These 
quantities are the first and second derivative of vryy and can be calculated 


according to (3.5). The trajectory position stsy is measured from t = 0, hence 
sty (t = 0) = 0. styy can be computed using (3.10) and (5.6). 


Unless stated otherwise, the set point of trajectory velocity vse and the set point 
of trajectory acceleration aset are defined as follows: 


Vset = VLim,TIY> ASet = 0 (5.11) 
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In contrast to (5.7), the optimization problem in (5.10) includes two additional 
summands that refer to the derivatives ayyy and jryy of the vrjy(t) trajectory. 
By penalizing deviations of arıy from aset and deviations of jryy from zero, 
the vryy(t) function is stabilized and uncomfortable driving is avoided, which 
can be caused by acceleration peaks and velocity oscillations. 


Each optimization goal has a corresponding weight. R7! denotes the weight of 
velocity error square, R7! the weight of acceleration error square and Rj! the 
weight of jerk error square. The reciprocals of the weights follow the interpre- 
tation of the filter algorithms from Subsection 3.3.2 and Subsection 3.4.2 by 
referring to the variances of the artificial measurements vset, dset and 0. R, is the 
variance of velocity measurement, R¿ the variance of acceleration measurement 
and R; the variance of jerk measurement. R} ! can be interpreted as a weight for 
the optimization goal of few travel time whereas both R7! and Rz! refer to driv- 


ing comfort. A suitable weighting combination reads R, = 5, Ra = 10, R; =1. 
This combination was derived by experiments, validated in real test drives and 
will be used in the remainder of this work. 


WLS and NWLS approximation algorithms usually solve an unconstrained 
problem and compute a function that is close to the data points. Problem 
(5.10) follows this approach and penalizes deviations of the function value 
from the data symmetrically. This interpretation of the data neglects that the 
data points of the upper speed limit also present a constraint to vryy. This 
constraint character is taken into account by using modified data in (5.10). Set 
point of trajectory velocity vser in (5.11) is not Vim Map, but the speed limit for 
trajectory optimization VLim Try: 


VLim,TJY = min (VLim, Map,’ [Pmin], + ++, VLim,Map,v [Pmax]) 
Pmin = round (styy — vryy ` AfLimayy + 1) (5.12) 
Pmax = round (styy + vtyy + AfLimTIY + 1) 
VLim,TJy İS the minimum of vLim.Map,» Within a spatial distance around the calcu- 


lated trajectory position styy. The distance depends on the trajectory velocity 
vryy and the temporal safety margin to upper speed limit Atrim try. AfLim.TIY 
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is a tuning parameter. In combination with At; sufficiently small it causes 
that there are enough data points at local minima of vLim,Map,v such that the 
trajectory does not exceed local minima of Vim Map,\- 


The nonlinear approximation problem (5.10) can be solved with LM and NRBA. 
To an approximation of (5.10) RBA can be applied. 


Advantage of solving the approximated optimization problem 


Regarding the control, [82] states that the time domain is beneficial because the 
vehicle dynamics itself is almost linear except from the air resistance. Therefore 
the problem can easily be simplified to a linear constrained quadratic program- 
ming (QP) problem. In contrast, the control in space domain is very nonlinear 
when the distance to the vehicle ahead is considered and therefore harder to 
be solved. Due to the spatially defined upper speed limit, (5.10) is a nonlinear 
problem though. 


However, RBA can solve an approximation of the nonlinear problem (5.10) by a 
quadratic problem. For this, in each iteration p the trajectory position stjy at the 
current time tp is determined by temporal integration of the trajectory velocity 
vryy with its currently estimated control point vector x in order to derive the 
upper speed limit vlim ty (Styy(tp)). Then vset = VLim,tyy is provided to the 
linear WLS algorithm RBA. The linear algorithm has no knowledge of the 
nonlinearity and assumes that the error vLim.rıy — Vryy depends linearly on x. 


The solution of RBA is not guaranteed to be optimal under these conditions. 
However, the effort is lower than with a method for NWLS and the solution 
of the approximated problem is satisfactory compared with an in general not 
globally optimal solution of the actual NWLS problem. For predictive vehicle 
control, [47, 28, 86] also approximated a nonlinear problem by a quadratic 
problem in order to be able to apply a quadratic programming method instead 
of a sequential quadratic programming method like LM. 


125 


5 Planning of velocity trajectories 


Influence of temporal safety margin to upper speed limit 


The upper diagram in Figure 5.2 depicts the velocity v versus the position s 
measured from the vehicle in driving direction. The dashed line shows the speed 
limit from map data with desired acceleration VLim,Map,v- The solid red line is 
a trajectory computed using RBA with Arıim.tıy = 18, At, = 1 s, Ath = 0.1 s 
and / = 1, whereby / denotes the number of spline intervals. Each red dot is 
a set point of trajectory velocity vse; that occurs during the iterations of RBA. 
For better view, only every tenth vse point is shown in Figure 5.2. The blue 
dots and lines show the same quantities for AtLim try = 4S. 


The temporal parametrization causes that in the v — s diagram the vse, data 
points are less close at higher speeds (e.g. for s < 200 m) compared to at 
low speeds (e.g. around s = 600 m). Due to the temporal safety margin to 
upper speed limit in (5.12), the spatial distance between the courses of Vset 
and VLim.Map.» is larger at higher speeds. Provided that Ary, is sufficiently small 
and that Arı;m,rıy is sufficiently large, there are enough data points to cause 
that the resulting trajectory does not exceed VLim,Map, at local minima, e.g. at 
s = 400 m and s = 600 m. 


Increasing AtLim,tyy also reduces the possibility of short-lasting velocity peaks. 
For example, at s = 500 m vLim,Map,; has a local maximum that both vset 
with Atlimtry = 1 s and the corresponding trajectory vrjy(t) reflect. For 
Attim try = 4 s, both vse and vrsy(f) cannot maintain the local maximum of 
VLim,Map,» at s = 500 m, however. 


Since the course of vse is not jerk-free and must be smoothed by the trajectory, 
close proximity of vryy(t) to vse: is not intended in all situations. The extent to 
which the trajectory can follow vse is determined by the weighting factors R;!, 
Rz; Rz! and the temporal distance of neighboring knots Ar,. Furthermore, the 
temporal distance of neighboring data points Aty needs to be sufficiently small. 
Reducing Aty increases the computational effort of RBA linearly. 


In the upper v — s diagram of Figure 5.2 vLim,Map,v is identical for all trajectories. 
The lower diagram depicts the same quantities versus time. In a v — ¢ diagram 
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For v - s diagram: --- VLimMap,» 
For AfLim,T1Y = 18: === VLim,Map,®  VSet 
For AfLim,T1Y =48: === VLim,Map,®  VSet 


Velocity in m/s 
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Velocity in m/s 


0 5 10. 15 20 25 30 35 40 
Time in s 


Figure 5.2: Influence of temporal safety margin to upper speed limit AfLim,Tjy on set point of 
trajectory velocity vse, and trajectory velocity vrsy. Only a subset vse, is shown. 
Upper diagram: The speed limit from map data with desired acceleration VLim,Map,v 1S 
identical for both trajectories. Lower diagram: vLim,Map,v differs with the trajectories. 
Adopted from [79]. 


VLim,Map,v Needs to be shown individually for each trajectory because when two 
trajectories differ, in general they have a different trajectory position at the same 
time t, hence vLim.Map,v also differs between the two trajectories. Adopted from 
[79]. A v- s diagram is often more convenient for comparisons because of the 
identical VL im Map,»- However, the lower v — t diagram in Figure 5.2 illustrates 
well that the temporal safety margin to upper speed limit leads to a constant 
time gap between the courses of vse and VLim,Map.v- 


Moreover, because of the constant temporal density of degrees of freedom 
resulting from the temporal distance of neighboring knots Atx, the curvature of 
the trajectory is independent of velocity in the v — ¢ diagram. In contrast, in the 
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— RBAp=100 —— RBAp=200 —— RBA,-300 
—— RBA ,=400 VLim,Map,v 


Velocity in m/s 


0 100 200 300 400 500 600 700 800 
Position in m 
Figure 5.3: Iterative trajectory optimization using RBA: Speed limit from map data with desired 


acceleration VL imMap,y (black) and trajectory after 100 (green), 200 (yellow), 300 (red) 
and 400 (blue) iterations p. 


v — s diagram the trajectory curvature usually decreases as velocity increases. 
Due to the equidistant knots with respect to time and (5.12), the trajectory 
can approximate the upper speed limit very close at low speeds because its 
degrees of freedom are dense with respect to position whereas at high speeds 
the approximation is less close. 


Advantage of iterative approach of RBA 


LM, RBA and NRBA compute the trajectory iteratively but in different ways. 
LM adapts all control points simultaneously, hence the whole trajectory is 
improved during an iteration. In contrast, RBA and NRBA process only a 
subset of the control points simultaneously and the temporal length of the 
trajectory increases with the number of iterations. 


Figure 5.3 shows a vrsy(f) trajectory generated using RBA after 100, 200, 300 
and 400 iterations with Ar, = 2 s, Atr = 0.1 s and Arıjm.rıy = 2 s. The iterative 
approach of RBA allows to pause the calculation and use the intermediate result 
before continuing the calculations. This possibility is beneficial in presence of 
restrictions on computation time. It is also known that the optimization is 


128 


5.3 Trajectory optimization 


finished for the subtrajectory that can be computed from knots and control 
points that are not used by RBA anymore but have been shifted out of the KF. 


Furthermore, with RBA the number of data points is not required to be bounded 
because of the shifting function. This is beneficial because a priori it is not 
known how much time data points tp the trajectory needs in order to reach the 
end of the electronic horizon. Also the length of the electronic horizon usually 
increases as the navigation system is recalculating the route (c.f. Section 5.1). 
RBA can start trajectory optimization immediately parallel to the route recalcu- 
lation and can add knots during the optimization process as needed. 


In contrast, LM cannot change the knots. Hence, the number of needed knots 
has to be overestimated, which causes additional computational effort, because 
more control point values are computed. 


Influence of temporal distance of neighboring knots and number 
of spline intervals 


Figure 5.4 shows trajectories with AtLim try = 1.5 s and Atr = 0.2 s. Trajecto- 
ries with temporal distance of neighboring knots At, = 1 s are depicted in red 
and trajectories with At, = 5 s in blue. For both cases trajectories computed 
by RBA with number of spline intervals J = 1 and J = 5 are compared to the 
LM solution. The solutions of both algorithms are similar and differences can 
mainly be identified from the acceleration diagram. 


Increasing the temporal distance of neighboring knots Ar, reduces the possible 
oscillation of velocity in a given time interval. This effect can be seen in the 
acceleration diagram for 200 m < s < 600 m. However, a large change of 
vry as for 600 m < s < 800 m can still be represented by At, = 5 s through 
correspondingly larger differences of temporally less dense neighboring control 
point values. For At, = 1 s, the RBA solutions are closer to the LM solution 
than for At, = 5 s. 
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VLim,Map,v 


--- RBAj-5 — LM 
--- BA — LM 


Velocity in m/s 


0 200 400 600 800 1,000 1,200 
Position in m 


Acceleration in m/s? 


0 200 400 600 800 1,000 1,200 
Position in m 


Figure 5.4: Influence of temporal distance of neighboring knots At, on trajectory optimization 
using RBA and LM and influence of number of spline intervals J on trajectory opti- 
mization using RBA. 


Increasing J also enables to smooth the trajectory because the filter can rework 
larger parts of the trajectory with hindsight when additional information is 
available. However, this can result in undesired velocity oscillations, e.g. for 
800 m < s < 1200s. This applies to both RBA and LM. 


Furthermore, RBA trajectories with larger J are more likely to exceed vLim,Map,v- 
This slightly occurs at s = 1200 m where the RBA trajectory with Ar, = 5s 
accelerates too early. The reason is a general flaw of the filter-based approach. 
When the filter adapts control points of previous spline intervals according to 
the knowledge of its state estimation covariance matrix, it does not compute any 
deviations or errors for these control points that could be penalized. Exceeding 
VLim,_Map,» can be avoided by increasing Atyim,tyy with Z. 
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As LM evaluates the whole sum of squared errors in the approximated NWLS 
problem in each iteration, this problem does not occur with LM. At around 
s = 1200 m both LM trajectories are identical. 


Although the above stated approach does not optimize the energy consumption 
by taking into account the electrical power explicitly, a reduction of energy 
consumption can be achieved implicitly by increasing the weight of acceleration 
error square R,'. The method of reducing acceleration and braking was also 
chosen in [1, 107, 112]. According to the investigations in [81] itis less effective 
than using an energy consumption model. Instead its advantage is that a simpler 
optimization algorithm can be used. 


5.4 Trajectory optimization considering 
electrical power 


O According to Section 4.4 low absolute values of the electrical traction power 


Pirac,elec lead to low power losses. Pirac elec is described by the adaptive electrical 
power model (AEPM) proposed in Subsection 4.7.2. The model output Pagem 
is given by 

Pagpm(t) = AEPM(vryy (t), ax). (5.13) 


The sensor longitudinal acceleration ax from (4.22) is computed as 
dx = arjy(t) + 8 - sin(a(stry(t))) (5.14) 


and the road slope angle «œ can be derived from the road slope stated in percent 
according to map data using (4.2). 


When considering the required electrical power in trajectory planning, (5.10) is 
augmented with the summand Ro. . (Pagpm)* which penalizes absolute values 
of PAEpM. Re denotes the weight of power error square and its reciprocal Rp 
is the variance of power measurement. 
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The resulting optimization problem reads 


P 


ĉ = arg min) (x, $ [VLim.TIY (stiy(tp)) = vovtp)] 
x p=1 


+R! > (dacs - arsy (tp)) + Rj! - (jrsv(tp)) 


+Rp' . (Pasow(tp))”) ; 


2 (5.15) 


As Pagpm correlates with the product of vryy and its derivative arıy, an approx- 
imation of this fourth optimization goal such that RBA can be applied, seems 
not promising for good results. Adopted from [79]. 


Most nonlinear filter algorithms such as UKF or PF do not accept a measure- 
ment matrix that relates a control point vector to a measurement as input. In- 
stead, they apply a given nonlinear measurement function to a set of control 
point vectors and use the resulting function values for estimating the relation 
between measurement and control points. However, for the first three optimiza- 
tion goals the exact values of the measurement matrix are known. These result 
from the values of the B-spline basis functions. 


Estimating the B-spline basis functions and additionally performing the tran- 
sition from a spatial to a temporal representation poses a challenge for the 
convergence of a nonlinear filter. Experiments were conducted with a square- 
root cubature Kalman filter (SCKF) stated in [6, 7]. The SCKF can be seen 
as a special case of the UKF. In experiments no weighting factor combination 
was found that yielded consistent convergence of the trajectory to the upper 
speed limit. This was also the case when the SCKF was applied to optimization 
problem (5.10), which does not consider electrical power. 


The MPF is beneficial for problems that can be subdivided in a linear subprob- 
lem and a nonlinear subproblem. A KF solves the linear subproblem optimally 
using a given measurement matrix while a PF is applied to the nonlinear prob- 
lem [155]. In the trajectory optimization application the nonlinear problem is 
the fourth optimization goal. The MPF converges for problem (5.15) reliably 
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provided that the weighting factors are within a certain range. For these reasons, 
NRBA includes a MPF. However, even though the MPF comprises a KF, the 
convergence radius of the MPF is still noticably smaller than that of a pure KF 
and the MPF is much more sensitive to filter parameter variations. 


O Braking for an upcoming curve usually requires negative Parpm for recu- 


. . . 2 
peration. Strong penalization of Pi ppm 


curves. In order to enable sufficient braking power, the computation of the error 


can prevent braking for upcoming 


e = VLimTJY — Vtyy Within the MPF is designed asymmetrically. If e < 0, the 
artificial measurement vLim try is replaced with vlimtry = Ry>vLim ` € + VIJY- 
Hence, e is multiplied with Ry>vLim if vrsy is above vlim try. In case of e > 0 
the original e is multiplied with Ry<vLim, the error weighting if vryy is below 
VLim,Tıy. In each case the resulting error square is afterwards weighted with 
R7! in the MPF. Adopted from [79]. 


Ry>vLim cannot be chosen arbitrarily large. First, Ry>yLim > 2 can prevent the 
MPF from converging to the upper speed limit. Second, during operation of the 
ALC situations can occur in which the vehicle exceeds VLim,Map.v. An example 
is that the driver intervens and accelerates in order to overtake a slower vehicle. 
At its beginning the trajectory needs to indicate the current vehicle velocity. 
Ry>vLim should be only so large that a clearly noticable braking occurs until the 
vehicle is below VLim,Map,v- For Ry>vLim > 2 braking according to the computed 
trajectory is uncomfortably abrupt. 


O The asymmetrical weighting of the velocity error can also be used in combi- 


nation with RBA but since there is no incentive to exceed VLim,Map,» because of 
the missing optimization goal regarding Parpm, the asymmetric weighting has 
very limited influence. Adopted from [79]. 


Influence of weight of power error square 


O By increasing the weight of power error square Rp or by lowering the vari- 


ance of power measurement Rp, respectively, the trajectory can be optimized 
towards low energy consumption. Figure 5.5 depicts various quantities for 
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—— NRBA, Rp = 10° —— NRBA, Rp = 10° —— NRBA, Rp = 10? 
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Figure 5.5: Relevant quantities of trajectories determined by NRBA and LM for different variance 
of power measurement. Velocity v, electrical traction power Pirac,elec, road slope 
y, energy loss between HV battery and wheel ELoss and HV battery energy Epa 


for trajectories determined by NRBA and LM that differ in the variance of power 
measurement Rp. Adopted from [79]. 
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trajectories determined by NRBA with Rp = 10°, Rp = 10° and Rp = 10%. 
For comparison, the LM solution for Rp = 10* is also shown. The remaining 
trajectory parameters are Af, = 2 s, Atr = 0.25 s, Afrimtiy =1sand/=1. 


In the depicted situation, all tractories indicate roughly the same velocity at 
the start position and also at the end position. Hence, the kinetic and potential 
energy of a vehicle at the beginning and end is roughly the same regardless of 
which trajectory it tracks. Therefore the energy consumptions resulting from 
following each of the trajectories can be compared. 


For Rp = 10* the depicted quantities for the NRBA and LM trajectories are 
very similar and both trajectories follow vLim,map,» Closely. With LM the energy 
loss between HV battery and wheel (ELoss) is ELoss = 77 Wh and the required 
HV battery energy (Egat) iS Egan = 721 Wh. For Rp = 10* NRBA achieves 
ELoss = 75 Wh and Egan = 714 Wh. 


When Rp is lowered, peaks in the electrical traction power Pirac,elec are reduced 
which translates to less energy loss between HV battery and wheel as well as 
less HV battery energy. Rp = 10° leads to ELoss = 71 Wh and Egan = 703 Wh. 
For Rp = 10? ELoss = 65 Wh and Egan = 684 Wh result. 


Regarding the driving style it can be observed that with lower Rp, the trajectory 
exhibits lower acceleration, stays more below VLim,map,», avoids short-lasting ve- 
locity peaks and tends to oscillate more. Adopted from [79]. The oscillations 
can mainly be explained by the larger slope variations which affect Ptrac elec but 
they partly also result from the fact that the MPF estimations are more strongly 
based on the state-space sampling of the PF than on the KF. With LM no en- 
ergy savings were achieved by lowering Rp, even after varying LM parameters. 
Instead, the trajectory diverged which indicates a narrow convergence radius of 
the unapproximated nonlinear optimization problem. 


O During the acceleration phase up to s = 1000 m the NRBA does clearly avoid 


peaks in Piracelec by reducing the acceleration at the slope. However, this does 
apply to the deceleration phase not to such an extent, even under consideration 
that from s = 2200 m on there is a high road grade which avoids the need for 
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large negative Pirac.elec. The filter could still smooth power demand more by 
decelerating earlier. It was also oberserved during road sections with negligible 
slope that NRBA decelerates too late and therefore the asymmetrical weighting 
of the velocity error was added. Adopted from [79]. 


As shown in Subsection 3.3.4 and Subsection 3.4.4, both novel filter based 
methods RBA and NRBA can reproduce a symmetrical approximation of WLS 
and LM, respectively. Provided that / is large enough, the filter lag is negligible 
with RBA. When the nonlinear optimization goal is weighted heavily with 
NRBA, a filter lag occurs, solutions become more oscillating and increasing I 
is not as unambigiously beneficial as with RBA. In the situation depicted by 
Figure 5.5, increasing J mainly leads to more oscillating trajectories whereas 
no significant change in deceleration can be noticed. Hence, the energy-saving 
potential of NRBA during recuperation is lower than expected. 


Increasing / should facilitate early deceleration because NRBA can adapt larger 
parts of the trajectory with hindsight when data that involves the lower speed 
limit is processed. MPF improvements that achieve the same approximation 
quality with less particles as stated in [203] would be beneficial for good results 
with larger J. They could help to overcome the flaw concerning deceleration 
because they enable to keep the state-space sampling density in each dimension 
constant when / is increased. 


At an early stage, at which the trajectory representation was not yet tempo- 
ral, the AEPM hyperparameters were not final and NRBA still included the 
SCKF instead of the MPF, automated tests for finding suitable target criteria 
weightings were done by A. Thorgeirsson [173] as a student assistant. 


5.5 Considering trajectory constraints 


After each iteration of the trajectory optimization, a lower limit on the trajectory 
velocity vryy is enforced and the trajectory acceleration arıy is restricted. This 
ensures that after any iteration a feasible trajectory is available. Due to the 
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high execution frequency the formulae are rather simple, avoid evaluating the 
B-spline function and take advantage of the convex hull property. 


vryy is restricted to nonnegative values, which translates to restricting all com- 
ponents of the state estimate X to nonnegative values: 


% = min (2,0) (5.16) 


Additionally the differences between neighboring components £ of X are lim- 
ited according to 


ĉi = max (ŝi + ddes,neg * At,, min En Xi1 + Ades,pos * At,)) > 


PS 2, 352 53 


(5.17) 


des,pos denotes the desired positive longitudinal acceleration and Ades.neg 18 
the desired negative longitudinal acceleration. The approach in (5.17) over- 
estimates the absolute value of trajectory acceleration using the convex hull 
property of the B-spline function. Hence, (5.17) might in some situations adjust 
control points although the specified acceleration limits are not yet exceeded. 


The following subsections state further constraints which are taken into account 
partly as hard constraints and partly as soft constraints. Hard constraints are 
constraints that reduce the solution space of an optimization problem and can 
lead to an empty set of feasible solutions. In contrast, the term soft constraints 
means that the optimization function is chosen such that violating these soft 
constraints leads to very bad function values with the aim that such solutions 
are avoided by the optimization method. 


5.5.1 Adaption to vehicle motion state 


O Small deviations between trajectory velocity vryy and vehicle velocity vyhel 


always occur because of an imperfect control. However, there are several cases 
in which both quantities differ significantly from each other. Examples are 
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driving off from a parking position and when the driver stops overriding the 
ALC by pressing accelerator or braking pedal. 


If the deviation exceeds a tolerated threshold, a new trajectory should be 
planned. The velocity and acceleration at the trajectory beginning should be 
adapted to the vehicle motion state given by vehicle velocity vype and vehicle 
longitudinal acceleration vvncı,x. Hence, the constraints 


vry (ti = 0) = Wnel (5.18) 
arıy(tı = 0) = Whel,x 


shall be enforced. Several methods for taking into account state constraints 
during Kalman filtering are presented in [163]. The state projection method can 
modify the control points such that the trajectory fulfills certain hard equality 
constraints at single points in time. Adopted from [79]. The state projection 
method projects an unconstrained estimate a onto the hyperplane defined by 
the constraints Dx = d. Usually Dx = d is an underdetermined equation 
system. The constrained estimate X, is given by 


-D7 (DD") (Dt; - d). (5.19) 


The difference of the projected control points need to be checked against the 
desired longitudinal acceleration. If they exceed the acceleration limits 4des neg 
OF Ades pos» the projection needs to be repeated with additional constraints that 
refer to differences between neighboring control point values. 


O The state projection method changes only the control points that influence 
the trajectory function at t; = 0. If the function is projected onto a vehicle at 
standstill and the temporal distance of neighboring knots Ar, between knots 
is small, the trajectory will demand an abrupt acceleration towards the upper 
speed limit unless control points that refer to following spline intervals are 
adapted too. 


A jerk-free and comfortable velocity transition can be achieved by modifying 
in the optimization problems (5.10) and (5.15) the desired acceleration value 
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Trajectory without adaption Trajectory with (1) 
—— Trajectory with (1) and (2) Trajectory with (1), (2) and (3) 
VLim,Map,v Vehicle motion state component 
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Figure 5.6: Individual effects of projection onto vehicle motion state using (5.19) (1), enforcement 
of acceleration constraints with (5.17) (2) and adaption of filter settings using (5.20) 
(3) on the trajectory during trajectory adaption to vehicle motion state. 


set (tp) from (5.11) and the acceleration weight R7! in the first filter iterations 
as follows: 


api Y. pi , 
Ra = Ran + Ra — Ran) /AtvnetAdapt * MiN(AtvhciAdapt 1), 


(5.20) 
Aset(tp) = Ývhncl,x + (0 — Vynci,x)/AfvhclAdapt > Min(AtvhelAdap 1). 


AtvnclAdapı and R3 A are tuning parameters. From an initially strong weighting 
Rh Of aset = Vvhelx att = O, the quantities R7! and aset return to their standard 
values in a linear fashion until to + AtvnelAdapt- Adopted from [79]. 


Figure 5.6 illustrates the individual effects of projection onto vehicle motion 
state with (5.19), enforcement of acceleration constraints using (5.17) and adap- 
tion of filter settings according to (5.20) on the trajectory. The trajectory param- 
eters are At, = 2 s, Atr = 0.1 s, Athimtry =2sand/=1. 
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In the depicted situation, the speed limit from map data with desired accelera- 
tion VLim,Map,» 18 constant at 26.7 m/s (black solid line in upper diagram). The 
vehicle motion state consists of the vehicle velocity vvncı = 20 m/s and the ve- 
hicle longitudinal acceleration Ývhci,x = —1 m/ s?. Each component is indicated 
by a black dot in the corresponding diagram. A trajectory without adaption (de- 
picted by green dotted line) starts at vLim,Map,» instead of vyncı. A trajectory that 
is only projected onto the vehicle state (orange dashed line) can demand strong 
accelerations since the projection only changes the first four control points. The 
red solid line results from additionally enforcing acceleration constraints with 
(5.17). The adaption of the filter setting using (5.20) with R; = 0.01 and 


a,to 
AtvnclAdapt = 150 s allows to get a more comfortable transition. 


The desired effect could also be achieved using a sigmoid function. However, 
the parameters of the sigmoid function need to be chosen carefully because 
filter sensitivity differs strongly with the magnitude of R,!. For simplicity a 
linear transition with rather high AfvhelAdapt Was used. 


5.5.2 Adaption to previous trajectory 


vry (tı), arıy(ti) and jryy(tı) with t; = O of a new trajectory can also be 
projected onto vryy, atyy and jryy of the current trajectory at its current eval- 
uation point f* in order to achieve a C? continuous connection between both 
trajectories. 


When connecting two trajectories using state projection, the knot positions 
of the new trajectory are additionally chosen such that the distance of tı to 
the neighboring knots equals the distances of t* to the neighboring knots of the 
previous trajectory. Then the shape of the new trajectory around 11 is identical to 
that of the previous trajectory around 1* and uncomfortable velocity oscillations 
at the joint are avoided. 
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Figure 5.7: Individual effects of projection (1) and improved control point initialization (2) on 


the new trajectory during adaption to current trajectory at evaluation point styy (t*) of 
current trajectory. 
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Figure 5.7 depicts various trajectories with Ar, = 1 s, AtLimtry = 1s, Atn = 
0.1 s and J = 3. At t* = 14.333 s or sqyy(t*) = 201.75 m, respectively, re- 
planning occurs and the new trajectory needs to be connected to the current 
trajectory depicted in green. The new trajectory depicted in orange is not 
projected and therefore the connection is not continuous. The new trajectory 
depicted in red is projected and a C? continuous connection to the current trajec- 
tory is achieved. However, shortly behind the connection point the courses of 
its derivatives deviate from the corresponding courses of the current trajectory. 
The blue trajectory results when additionally the control point vector of the new 
trajectory is initialized with the corresponding control point values of the cur- 
rent trajectory. The deviations in the courses of the derivatives are significantly 
reduced. 


However, in case of J = 1 the second derivative still deviates with this ap- 
proach. This is because the estimate cannot be improved with hindsight using 
the knowledge stored in the KF covariance matrix that results from additional 
data points. 


5.5.3 Adaption to vehicle ahead 


The Intelligent Driver Model (IDM) described in [174] provides an ACC func- 
tionality by computing a desired vehicle acceleration aıpm given by 


AIDM = Afree + Aint- (5.21) 


The free-road acceleration term afre With 


VSet 


ô 
v 
Afree = Ades,pos i z ES | (5.22) 


determines how the vehicle velocity vyhe converges to the set velocity Vset 
without a vehicle ahead. ages pos is the desired positive longitudinal acceleration. 
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The acceleration exponent 6 is usually 4. The interaction acceleration term dint 
with > 
Aint = Ades,neg (=) (5.23) 
AS VhclAhead 
causes the vehicle to decelerate in order to adjust the distance to the vehicle 
ahead ASvnclAhead to the desired minimum distance Asges given by 


Vvhel * (VVhel — VVhclAhead) 
Asges = ASJam + AtvnclAhead * VVhel + ` (5.24) 


2/Ades,pos i |@des.neg | 


ASVhclAhead and the velocity of the vehicle ahead vvnciAhead are Measurements 


of the radar sensor. Asjam denotes the desired traffic jam distance with a typical 
value of 2 m. AtvnciAhead is the desired time gap to the vehicle ahead and 
des,neg < O the desired negative longitudinal acceleration. 


Vset in (5.11) can be modified with aıpm during the iterative solution of (5.10) 
and (5.15) such that the trajectory optimization process takes into account 
the vehicle ahead as a soft constraint. The radar sensor data AsvhciAheaa and 
VvhclAhead refers to the new trajectory at start time t} = O. Under the assump- 
tion that vvnclAhead is constant and that the ego vehicle will follow the planned 
trajectory perfectly (vyna(t) = vryy (t), Vt > 0), AsvnciAhead can be calculated 
in following iterations as the sum of the initial distance AsypciAhead(fı) and 
VVhclAhead (t1) - t, from which the trajectory position sryy is subtracted: 


ASVhclAhead (tp) = ASvnelAhead(t1) + VVhclAhead(f1) * tp — STIY (tp) (5.25) 


In order to avoid dividing by small values, the IDM interaction free road term 
from (5.22) with velocity dependent ages,pos according to the acceleration tables 
and with vset = vryy is used only if vryy is larger than 2 m/s. For vryy < 2 m/s, 
free 1S set to the trajectory acceleration aryy. The interaction term remains as 
in (5.23). 


Vset is determined depending on a hysteresis with respect to din: If din > 
-0.05 m/s”, vser for trajectory optimization is computed according to (5.11). If 
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Figure 5.8: Effect of Intelligent Driver Model (IDM) on set point of trajectory velocity vse and 
resulting trajectory depending on the chosen desired time gap to the vehicle ahead 
AtvhclAhead during trajectory adaption to a vehicle ahead with a velocity of 12 m/s. 
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dint < -0.15 m/s?, braking because of a vehicle ahead is required and apm is 
used to correct vse, downwards: 


Vset = max (0, min (vpyy + Att, + dIDM; VLim,TIY)) (5.26) 


In experiments it turned out beneficial to let aset remain at zero. With these 
modifications the free-road term in (5.22) causes a comfortable acceleration or 
deceleration towards the current vryy coming from the upper speed limit. 


The minimum operation in (5.26) ensures that vryy can decrease fast enough. 
Without this minimum operation, the IDM can weaken the braking for a tight 
upcoming curve when there is no vehicle ahead. This is because (5.22) causes 
a deceleration towards a lower velocity set point vse that is sufficient and 
comfortable on highways but too weak for country roads, on which vse can 
vary strongly and deviations between vyhe and vset should be minimal. 


Figure 5.8 shows a situation in which the IDM adapts the trajectory using (5.26). 
If the vehicle trajectory tracks the planned trajectory and the velocity of the 
vehicle ahead remains constant, the specified desired time gap to the vehicle 
ahead will result. 


5.6 Scientific contribution 


DM are popular for automotive applications and have great potential in combi- 
nation with DP as their orthogonal features can complement each other. How- 
ever, the exponential growth of computational effort with increasing time hori- 
zon limits the application of DM to short time horizons. 


In general the resulting static optimization problem (2.7) in the DM approach 
is nonlinear and solved by sequential quadratic programming (SQP) techniques 
or interior point methods [187]. Frequently trajectory optimization is applied 
to vehicles with combustion engine or a hybrid power train, both of which have 
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Nonlinear problems, e.g. NWLS 


Convex problems, e.g. WLS 


Linear problems 


Figure 5.9: Classes and examples of optimization problems 


lots of degrees of freedom and constraints. For example, if the optimization 
problem incorporates the gear selection it is a mixed integer problem. 


In comparison, BEVs often have a power train with 1-speed gear box which 
simplifies the optimization problem. Quadratic problems are a subset of nonlin- 
ear problems that can be optimized more efficiently with QP methods because 
setting the derivative of the convex optimization function to zero results in 
sufficient conditions for the global optimum. 


The common linear WLS approximation problem is an unconstrained quadratic 
optimization problem and the NWLS approximation an unconstrained nonlin- 
ear optimization problem. Figure 5.9 depicts different classes of optimization 
problems. 


This work presents a DM based trajectory optimization approach for an ALC 
ofa BEV with fixed gear ratio such as the considered research vehicle. With 
the presented approach, the computational effort only grows linearly with the 
number of function parameters or the time horizon. This substantial saving is 
achieved by formulating the trajectory optimization problem as either aWLS 
approximation problem or aNWLS approximation problem. The approxima- 
tion problem is solved iteratively by RBA or NRBA, respectively. Each iterative 
method includes an iterative state estimator, also known as filter. 


Usually filters solve unconstrained linear and nonlinear WLS problems. There- 
fore limitations of the vehicle and restrictions of the environment enter the 
trajectory optimization problem in the presented approach mainly as soft con- 
straints. This means that the objective function is designed such that undesired 
solutions coincide with comparatively very bad values of the objective function. 
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However, some hard constraints are also enforced using the state projection 
method from [163]. Hard constraints narrow the solution space of the optimiza- 
tion problem. 


The trajectory is defined by a B-spline function and describes the desired vehicle 
velocity with respect to time. The trajectory results from an iterative solution 
of the approximation problem. Spatially or temporally defined soft constraints 
or optimization goals referring to derivatives of the function or its integral can 
be taken into account during the optimization process. For example, a spatial 
upper limit on the velocity defined by the trajectory results from the legal speed 
limit or tight curves. 


Treating such spatial constraints as soft constraints leads to aNWLS problem. 
However, the iterative trajectory optimization approach allows to approximate 
this problem in each time step by a WLS problem similar to [28, 47], who 
approximated a nonlinear problem by a quadratic problem in order to be able 
to apply a QP method instead of a SQP method. The linear state estimator for 
WLS problems enables large savings in computational effort compared to the 
nonlinear state estimator for NWLS problems. 


The energetic optimization of the trajectory is a case in which the approximation 
with a WLS problem cannot be expected to produce acceptable results. Then a 
nonlinear state estimator can be applied as demonstrated. 


The presented iterative trajectory optimization approach offers the following 
additional advantages: 


First, the number of function parameters does not need to be bounded. In- 
stead, the iterative estimator can determine additional function parameters if 
the temporal length of the trajectory needs to be increased. As this does not 
influence the computational effort in each iteration of the estimator, arbitrarily 
long trajectories can be planned. 


Second, the temporal length of the trajectory increases with the iterations and 
after each iteration an intermediate result is available for use. Therefore, the op- 
timization can be paused and continued in accordance with the time constraints. 
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Differentation from previous works 


The previous project e-generation [22, pp. 29-35] uses a spatially defined veloc- 
ity trajectory represented by a cubic polynomial. The polynomial is iteratively 
adapted to map data using a KF in such a way that an optimization problem simi- 
lar to (5.3) but with additional target criteria referring to the first two derivatives 
is solved. 


Due to few degrees of freedom of the cubic polynomial, only a trajectory with 
simple shape can be represented. The spatially defined velocity trajectory ex- 
hibits the flaws mentioned in Section 5.2. Especially the reluctant behavior at 
low velocities is noticable in test drives. A target criterion referring to the elec- 
trical traction power or the consumed energy is not included in the optimization 
problem. 


The approach presented in this chapter ensures nonnegative trajectory velocities 
at the spline function level by restricting its control point values via (5.16). In 
contrast, the nonnegativity of the polynomial is not ensured, only the result of 
the polynomial evaluation is restricted to nonnegative values. 


Furthermore, the previous approach adapts the trajectory only roughly to the 
vehicle velocity and acceleration by modification of the target values in the 
approximation, whereas in (5.18) an exact adaption via parameter projection is 
done. 


An adaption to a previous trajectory as in Subsection 5.5.2 is not needed with 
the polynomial function and constraints resulting from a vehicle ahead as in 
Subsection 5.5.3 were not considered. 
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Section 6.1 of this chapter states the stages of the development process of 
the ALC followed by a description of the architecture of the ALC and its 
components in Section 6.2. Section 6.3 investigates the influence of selected 
parameters on the ALC energy-saving potential in simulations. Section 6.4 
summarizes the technical contribution of this work concerning ALC. 


6.1 Development process 


The development of the ALC in the project e-volution followed the V model, 
which describes the software development process in the shape of the letter V 
and is depicted in Figure 6.1. The V model is widely used in the automotive 
industry [152, p. 25]. ISO 26262 states requirements for the development 
process of safety critical components and systems in vehicles. The procedure 
described in ISO 26262 is oriented towards the V model [187, pp. 110-111]. 


In the descending branch of the V model the customer requirements are ana- 
lyzed and translated into a logical architecture from which a technical architec- 
ture is derived. In subsequent steps this technical architecture is decomposed 
into systems and components. During this process on each level specifications 
and test cases are defined for later review of the development steps. The last 
step of the descending branch coincides with the first step of the ascending 
branch and includes the implementation of the specified components. 
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Figure 6.1: Software development process according to V model with allocation of in-the-Loop 
methods according to [187, p. 165]. 


The ascending branch comprises subsequent testing and integration steps from 
individual components up to testing the customer acceptance on the whole 
system. 


Each step of the descending branch defines specifications that have to be met 
during verification on the same level of the ascending branch as indicated by 
dashed arrows in Figure 6.1 [187, pp. 163-168]. 


The V model distinguishes four different test steps. The first three are compo- 
nent test, integration test and system test. These are verificiation tests. Verifi- 
cation denotes the process of evaluating whether an implementation meets the 
specified requirements for the corresponding development phase. Only the cus- 
tomer acceptance test is a validation test that determines whether all customer 
requirements are fulfilled. Methodological additions to the V model such as 
rapid prototyping enable validation at an early stage and help to avoid time- 
consuming and costly reworking loops across various stages of the V model 
[152, pp. 33, 152-154]. 
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Rapid prototyping denotes process steps that enable early verification or valida- 
tion of specifications. Typical process steps are modelling, simulation, integra- 
tion and test of the prototype in the vehicle [152, pp. 160-162]. Model-based 
methods can be applied for specification of software functions as well as vali- 
dation of specifications. Model-based software development allows to model 
software functions in a graphic environment using block diagrams [152, pp. 31- 
32]. A model compiler can compile the software functions for various target 
hardware such that the functions can be simulated as well as tested in the vehicle 
[152, p. 208]. If rapid prototyping methods support the development process, a 
virtual integration of the system is possible at the end of the descending branch. 
Before the integration steps this virtual prototype can be evalutated in simulated 
test drives [187, pp. 163-168]. 


In-the-Loop methods allow to combine models or real components at each de- 
velopment step with a reproduction of their real environment in order to obtain 
an assessable system. As there are no real components until the implementation 
phase of the V model, a simulation environment is used for virtual integra- 
tion. Figure 6.1 allocates different in-the-Loop methods to the steps of the 
development process. 


Model-in-the-Loop (MiL) enables a confirmation of the specification of cus- 
tomer requirements up to the logical architecture. The created model-based 
algorithms do not yet refer to the hardware of the target system. By transfer- 
ring the model-based algorithms into a simulation environment that is hardware 
independent but exhibits already technical characteristics similar to that of the 
target system one can perform Software-in-the-Loop (SiL) which allows for an 
assurance up to individual components. Hardware-in-the-Loop (HiL) denotes 
methods that transfer models from the SiL environment to real components. 
For example a function can be executed on a vehicle ECU while the system 
architecture of the vehicle is simulated. This allows to verify the interaction 
of the ECU with other simulated vehicle components [187, pp. 163-168]. HiL 
simulation usually requires real-time capable components [152, pp. 297-301]. 
When HiL methods are applied consequently, the entire system exists in real 
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components and can be verified up to the logical architecture. The Vehicle- 
in-the-Loop (ViL) method is especially useful for the development of driver 
assistance systems as it allows the operation of a real test vehicle in a virtual 
environment. Vehicle and virtual environment can be linked by either replacing 
the real sensors with interfaces to the virtual environment or by stimulating real 
sensors artificially. 


In the research project e-volution, the customer requirement was stated as an 
ALC that is as comfortable and reliable as its predecessor from the previous 
project and offers at least the same dynamic driving style as the ALC that 
was developed during the previous project e-generation for the same research 
vehicle. As main enhancements the e-volution ALC was required to be capable 
of planning of longer, more far-sighted trajectories and of realizing an energy- 
efficient driving style. 


The specifications of logical and technical architecture were taken over from 
e-generation. The logical architecture can be summarized as a system that 
controls motor torque according to map data. The technical architecture defines 
the interfaces of the ALC to the vehicle, for example which CAN messages are 
sent and received. Furthermore, it specifies the interfaces to the driver, e.g. how 
the driver can control the ALC using buttons next to the steering wheel and 
which information about the ALC state is displayed by the instrument cluster. 


System design specifies the software architecture that implements the ALC 
functionality. The functionality is divided into different components whose 
interfaces are defined. System design was conducted during e-generation and 
reused for e-volution. 0 The system design is specified by a framework 


in the model-based development environment MATLAB Simulink by Math- 
Works. The framework was created by an employee of the Porsche Engineering 
Services GmbH. It also takes care of communication with other vehicle com- 
ponents via several CAN busses and converts CAN messages into physical 
quantities. Furthermore, it checks the torque demand of the ALC system output 
for plausibility. This allows a safe vehicle operation in conjunction with an 
agile development of the ALC. Adopted from [79]. 
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Design of components converts the black-box description of the system de- 
sign into a white-box description for each module including the internal data 
flow such that each module can be implemented [187, pp. 168-173]. Within 
e-volution, component design led to the approximation algorithms RBA and 
NRBA and the AEPM. These components were implemented and tested in 
MATLAB. Then the trajectory planning based on RBA and NRBA was de- 
veloped and tested in MATLAB using the open-loop vehicle reference model 
mentioned in Section 4.5 to investigate its energy-saving capabilities. 


For further integration steps these functions were included into the Simulink 
framework. Translating MATLAB functions into sequences of Simulink blocks 
seemed not be advantageous. Therefore these functionalities were integrated 
into the Simulink framework as embedded MATLAB functions. As a result, 
each of the ALC components is based on few Simulink blocks that contain 
embedded MATLAB functions. Around the framework several simulink blocks 
including the closed-loop vehicle model from Section 4.5 were added in order 
to enable SiL tests of individual components and the whole ALC system. Ex- 
tensions for considering constraints stated in Section 5.5, such as the projection 
of the trajectory as well as taking into account the vehicle ahead, were almost 
exclusively developed within the Simulink framework because of the closed 
control loop. 


Additionally, a HiL test bench was used for component, system and integration 
tests as well as calibration. The HiL test bench was created by B. Fath during 
e-generation [22, pp. 60-73]. During HiL tests the developed functions run on 
the target hardware of the research vehicle. The target hardware is an ETAS 
ES910 rapid prototyping ECU and can execute compiled Simulink models. 
The HiL test bench allows to test the interaction of the software on ECU with 
the simulated vehicle architecture close to reality before the actual test in the 
vehicle. 


The HiL test bench consists of four hardware components. A real-time com- 
puter, a host computer, the prototyping ECU and a car PC. Figure 6.2 shows 
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Figure 6.2: Simplified illustration of HiL architecture according to [22, p. 66] 


the connections between the components of the HiL test bench in a simpli- 
fied manner. The real-time computer is an IPG Xpack4 that runs a simulation 
environment including all CAN communication between the ES910 and the 
architecture of the virtual vehicle. The host PC is an ordinary desktop com- 
puter that controls the real-time computer via an ethernet connection, starts 
simulations, records simulation data and provides a graphic simulation output. 
The simulation software is IPG CarMaker HiL, a vehicle dynamics simulation 
software that is extended by the capability to control the Xpack4. 


In CarMaker a vehicle model was created that is more detailed than the vehi- 
cle reference models described in Section 4.5 and also takes into account tire 
characteristics and vehicle kinematics. Furthermore, a model of the power train 
as described in Section 4.3 was integrated into the CarMaker vehicle model as 
a Simulink model and a road model of the WR was created from Open Street 
Map data and included in CarMaker. These models were also added by B. Fath. 


The host PC also runs the experiment environment ETAS INCA. Via a second 
ethernet connection of the host PC, INCA can flash compiled Simulink models 
onto the ES9 10, start them on the ES910 and view and change model parameters 
during run-time, e.g. for calibration purposes. 


The car PC is a small computer designed for automotive applications and con- 
nected to the ES910 via a CAN bus. It receives a GPS position that is generated 
by the simulation on the Xpack4 and passed trough the ES910. Using this 
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information the map software EB Assist ADTF by Electrobit, which runs on 
the car PC, provides map data for the ES910. 


A similar HiL setup also for a trajectory optimization and vehicle control use 
case is described in [54, pp. 181-192]. 


The HiL environment enables safe, reproducible and automatically performed 
tests. It also helps to reduce expensive testing time in the real vehicle because 
most specifications can be verified in advance [187, pp. 163-168]. 


Despite the great benefits of the in-the-Loop methods, they cannot completely 
replace real test drives. At further development stages new test scenarios that 
cause problems more likely occur in reality than in the already well-known 
and simplified simulation environment. Furthermore, some features of driver 
assistance systems require a subjective assessment [187, pp. 163-168]. 


O Real test drives were conducted on the WR itself as well as on the highway 


and at roundabouts nearby. The test drives could focus on determining various 
calibration parameters such as acceleration look-up tables, weighting factors of 
the trajectory optimization and parameters of the controller that computes the 
torque demand. Parameter settings for three different modes that correspond to 
different driving styles were derived. In a final acceptance test about ten em- 
ployees of the project partner who work in automotive research and engineering 
tested the ALC on self chosen routes. Adopted from [79]. 


6.2 System architecture and design 


O The software components of the ALC are denoted modules. The ALC system 


consists of route data module, parameter adaption module, trajectory module 
and controller module. Figure 6.3 depicts the architecture of ALC. The fol- 
lowing sections describe the component design, i.e. the function each module 
performs. Adopted from [79]. 
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Figure 6.3: System architecture of the automated longitudinal control. Adopted from [79]. 


6.2.1 Route data module 


O The input to the route data module comprises six vectors that contain the 
values and positions of legal speed limits, the road curvature and the slope for 
the road section ahead of the vehicle in driving direction. Similar to a run- 
length encoding, the vectors contain only information for positions at which the 
corresponding quantity changes significantly. The map data is provided via the 
CAN bus by the navigation system. Adopted from [79]. Since no destination 
is specified, the navigation system only knows the route up to the next junction 
or roundabout. If all roads connected to the junction have the same average 
amount of traffic, the electronic horizon ends at this junction. If one road has a 
higher traffic density, the navigation system assumes that the driver takes this 
road because it is the most probable path (MPP). 


Following the MPP repeatedly leads onto the next larger road and finally onto a 
highway. Unless the driver follows an assumed MPP, there will be a short time 
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period with no available map data. The navigation system then recalculates a 
new path whereby the length of the electronic horizon quickly increases. 


Sometimes the driver likes to turn at a junction but the navigation system as- 
sumes that the junction will be crossed straight. In such cases the ALC does 
not reduce the vehicle velocity sufficiently and the driver needs to brake. Since 
the ALC neither uses information about right of way nor can detect traffic signs 
or traffic lights via a camera, the driver usually has to override the ALC at 
junctions and roundabouts. 


The length of the electronic horizon ranges from several hundred meters in 
urban areas over roughly | km on country roads to several kilometers on high- 
ways. 


The route data module is executed every 500 ms and mainly converts the map 
data to a meter discrete representation for up to 3 km ahead of the vehicle 
position and provides three vectors that describe the courses of road slope, road 
curvature and legal speed limit. However, the input vectors of the legal speed 
limit can be modified before the conversion by three mechanisms: 


First, the driver can overwrite the current legal speed limit of the map data by 
adding a positive or negative offset using buttons next to the steering wheel. 
This adaption functionality proves beneficial in situations in which the legal 
speed limit in the map data does not match the actual legal speed limit, for 
example at temporary construction sites. 


Second, the positions of upcoming speed limits can be changed. If the legal 
speed limit is increased, the ALC will not start to accelerate the vehicle before 
the new traffic sign is passed because the planned trajectory does not exceed 
the generated upper speed limit. When leaving a city, this behavior is usually 
undesired by the driver and can provoke overtaking maneuvers of following 
vehicles. Therefore two parameters were added. Both describe a time offset, 
one for the case of an increase of the legal speed limit and one for a decrease of 
the legal speed limit. In both cases the corresponding parameter is multiplied 
with the lower of the two speed limits to determine the spatial distance by which 
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the position of the higher speed limit is shifted in the map data. This allows to 
implement a more realistic driving behavior or adaption to driver preferences. 
For example, it can be achieved that the ALC already starts to accelerate when 
the driver sees the traffic sign that indicates a higher legal speed limit. 


Third, an upper limit is enforced on the legal speed limit depending on the 
driving mode of the ALC selected by the driver. This is needed on highways 
because there is not always a legal speed limit. 


The route data module was implemented by D. Dörr [22, pp. 23, 24]. The 
author of this work added the time offset parameters and logic for processing 
roundabout data. 


6.2.2 Parameter adaption module 


O The parameter adaption module updates the ATFM and AEPM described 
in Subsection 4.7.1 and Subsection 4.7.2, respectively, every 50 ms using ve- 
hicle data from the CAN bus. During the ATFM update one KF iteration is 
performed and during the AEPM update the FB-KRLS takes into account a new 
measurement. 


An update of the ATFM requires the traction force Firac, the vehicle velocity 
Vvhe and the sensor longitudinal acceleration a, and an update of the AEPM 
the electrical traction power Pirac.elec» VVhcl and ax. While vyncı and a, are mea- 
surement signals on the CAN bus, Frac needs to be approximately computed 
from the motor torque signals using (4.25). Pirac,eiec can be determined from 
the voltages and currents of the electric motors that are available on the CAN 
bus. 


Before CAN signals are fed to an estimator or used to calculate the input quan- 
tity for an estimator, signal noise is removed using a polynomial function ap- 
proximation method [145], also denoted polynomial Kalman smoother [144]. 
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In two cases the models are not adapted: First, during a vehicle standstill 
because the missing excitation in the data can cause estimations to diverge 
from their true value. Second, if CAN signals indicate that there is signifi- 
cant hydraulic brake pressure because the models only consider braking via 
recuperation and are not valid when hydraulic brakes are active. 


Output of the parameter adaption module are updated parameters of ATFM and 
AEPM. The parameter adaption module including ATFM is from F. Bleimund 
[22, pp. 24-29]. The author of this work added the AEPM. Adopted from [79]. 


6.2.3 Trajectory module 


O The trajectory module performs the tasks that Chapter 5 described in detail. 


These tasks include generating an upper speed limit and planning a trajec- 
tory. Inputs are the map data from the route data module, parameters of the 
AEPM provided by the parameter adaption module, vehicle data and the se- 
lected driving mode. The selected driving mode determines which longitudinal 
acceleration look-up table is used to compute the speed limit from map data 
with desired acceleration. Furthermore, in driving mode "Normal" and driving 
mode "Sport" RBA performs the trajectory optimization whereas in driving 
mode "Range" NRBA optimizes the trajectory using the AEPM, whereby the 
absolute value of the electrical traction power is penalized strongly. 


The trajectory module is called every 500 ms and then can either continue 
planning the current trajectory or start planning a new trajectory. 


Continuing planning the current trajectory means that the module performs 
additional iterations of RBA or NRBA until it reaches the iteration limit for the 
current module call or the optimized trajectory reaches the end of the electronic 
horizon (c.f. Figure 5.3). 


Planning of a new trajectory is only triggered if distance limit or time limit since 
the last beginning of a trajectory planning have been exceeded, if the driver has 
stopped overriding the ALC by pressing the braking pedal or accelerator pedal 
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or if data of the curvature vector has changed more than a threshold since the 
last module call indicating that the driver has chosen a route that differs from the 
MPP. Additionally, the length of the electronic horizon must exceed a specified 
lower limit to enable planning of a new trajectory. Adopted from [79]. 


In two cases, vryy and aryy of the new trajectory are projected onto the vehicle 
velocity vyhe and vehicle longitudinal acceleration Vye] x according to Subsec- 
tion 5.5.1: First, if there is no valid previous trajectory, e.g. because map data 
was temporarily not available, and second, if the current vehicle velocity devi- 
ates more than a certain threshold from the planned velocity vryy, e.g. because 
the vehicle is starting from a standstill or because of an intervention by the 
driver. 


In all other cases vryy (t1), arıy(tı) and jryy (41) with t; = 0 of the new trajec- 
tory are projected onto vrsy(f*), arıy(t*) and jryy(t*) of the current trajectory 
at its last evaluation point t* according to Subsection 5.5.2 in order to achieve a 
smooth connection between both trajectories at £*. 


The knot vector k and estimated control point vector X are the outputs of the 
trajectory module. 


6.2.4 Controller module 


O Every 20 ms the controller modules computes the desired velocity Vaes and 
desired acceleration ages and translates these quantities into a motor torque 
demand Ties that causes the vehicle to track the velocity trajectory. In most cases 
Vdes and Ages equal the trajectory velocity vryy and the trajectory acceleration 
arjy, respectively. vrsy and arıy are determined by evaluating the trajectory 
defined by the knot vector and estimated control point vector from the trajectory 
module at the current point in time measured starting from the last trajectory 
planning start. Adopted from [79]. 


However, the controller also comprises three functionalities of the trajectory 
module: 
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1. Trajectory adaption to vehicle motion state: While the driver is overriding 
the ALC by using the pedals, the controller module projects the trajectory 
onto the vehicle motion state given by vehicle velocity and vehicle accel- 
eration as stated in Subsection 5.5.1. Therefore, as soon as the ALC takes 
over control again, the controller can immediately evaluate the trajectory 
that is connected jerk-free to the vehicle motion state. During the next 
call of the trajectory module planning of a new trajectory is initiated. 


2. Trajectory adaption to the vehicle ahead (c.f. Subsection 5.5.3): The 
IDM in the trajectory module serves mainly for planning the approach 
to the vehicle ahead. For the actual vehicle following task an IDM is 
implemented in the controller as well. In order to be able to modify vryy 
and arıy according to IDM, the controller module can perform the state 
projection method described in Subsection 5.5.1. 


3. Enforcement of velocity and acceleration constraints: If the controller 
changes the control points of the trajectory, it must also enforce the 
restrictions vryy > O by adapting control point values and enforce re- 
strictions with respect to aryy by adapting the differences of neighboring 
control point values (c.f. Section 5.5). 


Since the controller module can alter the trajectory, the trajectory module ac- 
cepts a control point vector coming from the controller module as input and 
interprets it as the previous trajectory. If there is no valid trajectory available 
for evaluation, the controller module causes the desired vehicle acceleration to 
quickly diminish and keeps the velocity constant. If no map data is available 
for several seconds, the ALC system deactivates itself. 


The in this subsection above stated functionality was added by the author of this 
work. However, the architecture and control loop including the pilot control 
were adopted from F. Bleimund [22, pp. 35, 36]. 


O Figure 6.4 depicts the controller module architecture and the control loop. 
Vdes and Ages are inputs to a pilot control and a model predictive control (MPC). 


161 


6 Automated energy-efficient longitudinal control 


Pilot control Tdes,PC 
based on ATFM 


Vdes; Ades 


VVhch VVhel,x 


Figure 6.4: Architecture of controller module and control loop. The pilot control generates an open- 
loop motor torque demand Tyes pc based on desired velocity Vdes, desired acceleration 
des and road slope angle « using the adaptive traction force model (ATFM). The 
model predictive control (MPC) computes a closed-loop torque demand Tyes mpc to 
minimize the remaining deviation of vehicle velocity vyhe and vehicle longitudinal 
acceleration Vyhe]x from Vges and Ages. Adopted from [79]. 


The pilot control contains the ATFM and computes an open-loop torque de- 
mand Tdes pc based on road slope angle œ, Vaes and dges. Due to imperfect 
map data, sensor data and ATFM, the vehicle velocity vyn.ı deviates from vdes 
and the vehicle longitudinal acceleration vvncı,x deviates from dges. The two- 
dimensional MPC computes a closed-loop torque demand Ties mpc in order to 
minimize these deviations. 


The torque demand of the controller module Tyes is the sum of the pilot control 
torque demand 7ges,pc and the MPC torque demand Tyes mpc and reaches the 
motor ECU via the CAN bus. The motor ECU computes the distribution of the 
torque demand between the front and rear motor using the strategy described 
in Section 4.3 and controls these actuators accordingly. If the torque demands 
indicates a deceleration request, it can additionally be allocated to the brake 
ECU that controls the hydraulic brakes. 


The control loop is closed through the vehicle where CAN messages transfer the 
feedback information. A comprehensive review of trajectory tracking methods 
is provided in [40]. 
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Pilot control 


The pilot control determines an open-loop torque demand Taes,pc With 


Frac Ydyn, FA + Fdyn,RA 


Taes,PC = (6. 1) 


ig 2 

so that the vehicle roughly tracks vaes and Ages. ig is the gear ratio, Fayn,ra the 
dynamic front wheel radius and rayn,ra the dynamic rear wheel radius. The 
traction force Frac results from the ATFM in Subsection 4.7.1: 


Frac = (1, ax, Vies) "X Vhcl (6.2) 
— 


=:C vnel 


The parameter adaption module provides the adapted vehicle parameter vector 
Xvhel- The vehicle motion vector Cvncı depends on Vdes? and ax, whereby 


Ax = Ades + g Sin(a@). (6.3) 


Adopted from [79]. According to Subsection 4.2.1, the road slope angle œ can 
be computed from the road slope y, which can be derived from map data or the 
signal of a slope estimator that is available on the CAN bus. 


The simple calculations in the pilot control are very transparent, the estimated 
vehicle parameters can easily be limited to a valid range and the other values 
are measurements. The pilot control contributes to the torque demand to a large 
extent and serves for overcoming the majority of the sum of driving resistances. 
Using the pilot control allows to restrict the output of the MPC that closes the 
control loop to low absolute values without risiking a vehicle standstill at a 
slope. Restricting the MPC output seems beneficial from safety considerations 
because the MPC performs operations that can be numerically problematic such 
as dividing by potentially small values. Additionally it allows to apply a simpler 
vehicle model in the MPC for less computational effort. 
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Model predictive control 


Model predictive control (MPC) is a direct trajectory optimization method 
that approximates the constrained infinite dimensional OCP from (2.1) by a 
finite dimensional problem that can be solved efficiently in real-time. The 
approximation occurs by discretizing the problem and solving it for a finite 
time horizon T. 


MPC uses a system model to predict the future system behavior. A possible 
linear time invariant system model without uncertainties reads 


Xp+1 = Ax, + Buy (64 
Yp =CXp 


with system state xp, control input up, system output y, and discrete time 
index p. In each time step t, MPC solves the following optimization problem 
for the horizon t = 0, öt,...,T with cycle time ôt: 


a, = argmin J (x,,u,) 
Ur 


Txt Us) = 167 - Cxi) Q (y: - Cxi) +u, Run] (65) 
i=0 


subject to Fx, + Gu, <1 


Herein the control input sequence u; = {uojr, Uilt, .. ., UT—1]1) Over the horizon 
is the optimization variable. x;¡, and w;j- denote the predicted values of x and 
u, respectively, for time f + i based on the information available at t, whereby 
Xojr = x; is assumed. Present and future set points are specified via the system 


output y; = [Yojr Y 1]r ---> Y T-111). 


The cost function J includes symmetric error weighting matrices R and Q, 
which are assumed positive definite and positive semidefinite, respectively. 
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ot time horizon T 
A | > 
past | future 
> 
A t 
NL 
úl 
| | | 
| | | gi 
fp-3 tp-2 tp-1 tp t 


Figure 6.5: Model predictive control optimization horizon with time horizon T, cycle time ôt, 
current time step tp, planned and actual state X, x, planned and actual control input @, 
u [187, p. 1431, adapted]. 


The matrices F and G allow for linear constraints with respect to x and u, 
respectively. The constraints are formulated as inequalities, which apply ele- 
mentwise. 1 is a matrix of ones. For increased robustness and in order to avoid 
an empty set of feasible solutions, constraint softening can be applied. Thereby 
hard constraints become penality terms in the cost function analogously to the 
approach stated in Section 5.5 [61]. 


Of the optimal control sequence û, only the first element #9); is applied to the 
system. At the following time step the process is repeated with the additional 
information x;+] on the system state. 


Figure 6.5 illustrates this MPC approach, that creates a feedback which partially 
compensates for model inaccuracies and leads to a closed-loop control that can 
take anticipatory control actions for future set points and events. The abilities 
of looking ahead and considering constraints are also distinguishing features of 
MPC compared to a PID controller [187, pp. 1431-1432], [95, pp. 2, 13-16]. 
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An analysis of different controller designs is provided by [169] and different 
vehicle models for controller design are stated in [136, p. 12]. 


Vehicle related MPC applications in literature include ACC systems, e.g. [13, 
61, 112, 154]. These use MPC as a high level controller that calculates a 
sequence of acceleration commands for the ACC equipped vehicle so that it 
keeps the required distance to a vehicle ahead, whereas a low level controller 
tracks the acceleration sequence using commands to power train and brakes. In 
the MPC the vehicle is represented with a first-order model because the low 
level controller cannot realize commands instantaneously due to the limited 
bandwidth of the vehicle. 


In this work the MPC serves only as a low level controller and is additionally 
supported by the pilot control that considers the driving resistances. Further- 
more it seems reasonable to assume that the research vehicle is capable of 
following power train and brake commands relatively quickly because of the 
high power to weight ratio and the missing delay from changing gears. For 
these reasons a vehicle model similar to the one in [116] is used. It is not a 
first-order model but a double integrator involving vehicle velocity vype, vehi- 
cle longitudinal acceleration vyncı,x and vehicle longitudinal jerk Vvncıx. The 
forward difference approximations 


Whel (t + ôt) — Vynel(t) 


Vvholx(£) = 

SER (6.6) 
m Vvhel,x (t + ôt) - Vvhel.x (t) 
Ÿvhcl,x (t) = St 


give the time-discrete state-space representation of the vehicle model, in which 
VVhcl,x is the control input [13, 116]: 


VVhel,1+1 1 6t\ | Vvnel,: ig 0 E (6.7) 
= : * VVhel,x,t . 
VVhclx,1+1 0 1 VVhel.x,t 1 
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The system state comprises vehicle velocity and acceleration and should be 
equal to the set points for the desired velocity vge, and desired acceleration ages. 
Hence, the quantities in the state-space model of (6.4) read: 


1 ot 0 1.0 
A = 5 B = 7 C = 5 
0 1 1 0 1 
a VVhel Vdes 
Ut = VVhelx» Xt =]. > y= 
VVhcl,x Ades 


Simplifying Firac in (6.1) to Frac = Myhel * VVhcl,„ and calculating the temporal 


(6.8) 


derivative shows that the control input u defined as the vehicle longitudinal jerk 
Vvhelx translates to changes of MPC torque demand, Ties mpc: 
OTaes MPC  Mvyhcl * Vvnelx Fayn,FA + Fdyn,RA 
FF Or (6.9) 
dt 1G 2 
Discretization of the above equation using the backward difference approxima- 
tion gives Taes. mpc for the current time step: 


Myhel * Of Yayn,FA + Fayn,RA 


io 5 (6.10) 


Taes,MPC,t = Taes,MPC,t-1 + VVhel,x,t * 


Thereby Ties, mpc,o = 0 and ôt = 0.2 s are used. The change of Tues mpc in each 
ôt is limited to [AT min, AT max] = [-50 Nm, +50 Nm] with following matrices 


in (6.5): 
0 0 = Myhel -Ot-(Tayn,FA FFdyn,RA) 
= = 2-ig AT min 
Fe 00 G= Mvyhcl Ôt: (Tayn, FA +rdyn, RA) (6.11) 
2-iG AT max 


Furthermore, the R describing the cost of the jerk and Q describing the cost 
of deviating from the trajectory velocity and the trajectory acceleration on its 
main diagonal elements are set to 


R=10%, Q = (20 1)-1. (6.12) 
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The horizon comprises the current time step £ and nine future time steps, which 
results in 1.8 s look ahead of the MPC. With a previously in e-generation 
applied PID controller uncomfortable oscillations between acceleration and 
deceleration occured on test drives while driving downhill with roughly constant 
velocity. After the PID controller had been replaced with the MPC, these 
oscillations were not noticed anymore. 


However, the computational effort of MPC is usually much larger than that 
of the PID controller and depends on the system model and the prediction 
horizon [147]. Problem (6.5) has a convex quadratic objective function and 
linear constraints, which allows using various QP solvers based on active set 
methods or interior point methods for performing the online MPC optimization. 


This work uses the general purpose QP solver provided by [94] with minor 
modifications that consist of rewriting instructions that are incompatible with 
compilation for the target hardware and of removing functionality and depen- 
dencies that are not needed to solve the above stated problem. Modifications 
of the solver, first implementations and tests of promising state-space models 
within a provided simplified control loop were done by C. Lee [98] during his 
time as student assistant that started after completion of his Master’s thesis. 
Parameter tuning and closed-loop testing within the ALC in simulated and real 
test drives were done by the author of this work. 


The disadvantage of general purpose QP solvers is that they do not exploit the 
special MPC problem structure and therefore using them might prevent the ap- 
plication of MPC to problems with a high sample rate, high-dimensional model, 
or long time horizon T. For example, the computational effort of both active 
set and interior point methods increases roughly cubically with T, whereas the 
effort of algorithms that exploit sparse matrices increases only linearly [95, 
p. 42], [147]. A survey of technologies for linear and nonlinear MPC is given 
by [137] and [171] provides a comparison between PID controller, linear MPC 
and nonlinear MPC. 
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6.3 Energy-saving potential and effects of 
parameters 


O This section investigates the energy-saving potential using the automated 


longitudinal control (ALC) with different parameter settings compared to drives 
with manual longitudinal control (MLC) on the Weissach route (WR) defined 
in Section 4.6 in simulations. 


Automated longitudinal control setup and data 
generation 


The varied parameters are the temporal distance of neighboring knots Ar,, the 
number of spline intervals /, the temporal safety margin to upper speed limit 
At jm,rıy and the weight of power error square Re 


Simulations of the ALC are conducted for all possible parameter combinations 
that result from At, = {1.5 s,3 s}, Z = {1,3}, AtLimtry = {1 s,2 s} and Re = 


1 1 i At. O! 
{ 10000° 5000° 1000» 500° Too 50}. Adopted from [79]. 


During a simulation of the ALC on the WR, the position s on the WR is calcu- 
lated by integrating the vehicle velocity over time. With s the corresponding 
road curvature and road grade are determined from the map data of the WR. 


Section 5.4 applies the open-loop reference model from Section 4.5 in a back- 
ward simulation in order to calculate the energy consumption that results when 
the vehicle tracks a planned velocity trajectory perfectly. The current section 
considers the whole ALC including its controller module and calculates the 
energy consumption that results from the actually realized velocity trajectory 
using the closed-loop reference model from Section 4.5 in a forward simulation. 


O In order to compare different ALC settings, the required trip time fyyip for 


completing the WR is converted into the average velocity by dividing tTrip 
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ALC with Atpim try =2: ~O- Ate = 1.5,1=1 --@- Ate =31=1 

ALC with Atyimtyy = 2: ~= At, =15,1=3 a Atk = 3,7 = 3 

ALC with Atuim try = 1: -G- At, =15,1=1 -@- At =3,1=1 

ALC with Atuim try = 1: -0- At, =1.5,2=3 -4- At, = 3,7 =3 
MLC: —— 
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Figure 6.6: Average energy consumption vs. average velocity on the Weissach route with auto- 
mated longitudinal control (ALC) under different parameter settings in comparison to 
adapted and resimulated real drive with manual longitudinal control (MLC). Adopted 
from [79]. 


by the length of the WR and the energy consumption is scaled to the energy 
consumption per 100 km. This approach allows to summarize each simulation 
in a single data point in Figure 6.6, which depicts energy consumption versus 
average velocity. The lines are approximations of data points differing only in 
RP by quadratic polynomials. Adopted from [79]. 


For some parameter combinations the values of RE are high or low enough to 
make the MPF instable. In such cases the resulting energy consumption is very 
high while the average velocity is either very low or very high. Such data points 
are excluded from Figure 6.6. Reduced convergence radii are noticed mainly 
for the combinations that include Arıim.rıy = 2 s and At, = 1.5 s. AtLim,TJY 
influences the data provided to the MPF and for nonlinear filters the covariance 
also depends on the data. If Ar, is reduced while Aty is kept constant, the 
covariance matrix entries tend to remain larger. This suggests that for these 
settings, the MPF parameter Q4, Q need to be reduced. 
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Manual longitudinal control setup and data generation 


O For comparison with a drive with MLC, a real test drive on the WR with 
MLC was performed and CAN data was recorded. From the recorded data, only 
the vehicle velocity vyncı with a sample rate of 100 Hz was further used. By 
integrating vvncı Over time, the corresponding position on the WR for each data 
point is computed. Time information is then discarded because the following 
modifications of vyncı do not depend on time: 


First vyncı is multiplied with a factor around one in order to achieve MLC drives 
with various average velocities in the same range as the average velocities of 
the ALC simulations. 


During the real test drive the driver had to stop the vehicle at junctions in order 
to check for other traffic participants with right of way. When the vehicle 
has to come to a standstill and drive off again this increases both the energy 
consumption and time required for completing the WR. 


The simulation environment neglects right of way and other traffic participants 
do not occur. Therefore the vehicle passes junctions with a low speed without 
coming to a standstill. In order to not disadvantage MLC drives in the compari- 
son, position dependent lower limits on vyncı are applied to MLC drives. These 
lower limits coincide with the velocities the ALC chooses at the corresponding 
junction. Varying ALC parameter settings influences the velocity with ALC at 
the occuring junctions only slightly. 


During the real test drive there was few traffic and only once during a short 
period the driver was restricted in the choice of vehicle velocity because of 
a vehicle ahead. At the corresponding position the recorded vyhc data was 
adapted to a typical ALC behavior. 


The ALC computes an upper speed limit for the route ahead, which among 
others depends on the specified maximum absolute value of lateral acceleration 
Vvnely,max- The originally recorded vyna results in a lateral acceleration larger 
than the maximum absolute value of lateral acceleration for the ALC. Increasing 
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WVhcl,y,max leads to higher average velocity and less energy consumption. In 
order to compare ALC and real drive under same conditions, an upper limit 
is enforced on the vyhe of the MLC drives. The upper limit is the speed limit 
from map data with desired acceleration. 


After multiplying the recorded vype1 with a factor around one and imposing the 


various limits the corresponding time £ is re-calculated using t = = F and the 
MLC energy consumption is determined with the open-loop reference model. 
Figure 6.6 depicts the resulting data points and quadratic approximation for the 


MIC in black and by a solid line. 


The quadratic approximation line serves as a benchmark. The ALC data points 
should be below this line, meaning that the ALC achieves the same average 
velocity with less energy consumption or a higher average velocity for the same 
energy consumption. 


Initial evaluations of the energy consumption of MLC and ALC in a non-final 
setup were performed by A. Thorgeirsson [173] as a student assistant. 


Parameter effects 
Temporal distance of neighboring knots 


Approximations for data points that differ only in the temporal distance of 
neighboring knots Ar, and weight of power error square RE indicate higher 
average velocities and mostly also higher energy consumptions for Ar, = 1.5 s 
than for At, = 3 s. The reason for this is that trajectories with larger Ar, cannot 
follow the upper speed limit as closely because their degrees of freedom are 
temporally less dense. Hence, the vehicle is generally slower and omits some 
inefficient velocity peaks. Adopted from [79]. 


172 


6.3 Energy-saving potential and effects of parameters 


Number of spline intervals 


When the number of spline intervals / is increased, the MPF can adjust the 
trajectory in more previous trajectory intervals. When Re is chosen high for 
strong penalization of the electrical traction power, a large J should be beneficial 
because the MPF can for instance retrospectively reduce the acceleration if the 
road grade suddenly increases and maintaining the acceleration would lead to 
large power demand and therefore large power loss. 


However, similar to Chapter 3 and Chapter 5 the simulation results do not 
support this statement. All approximations for Ar, = 1.5 s and J = 1 are below 
those for At, = 1.5 sand J = 3. For At, = 3 s the ranges of the average 
velocities of the approximations for J = | and J = 3 do not always overlap 
but they can still be compared using their positions with respect to the MLC 
approximation. Using this criterion the approximations for J = | are always 
better. 


The reason for the worse results with larger J can be explained by the fact that 
the PF in the MPF uses a set of particles, each of which is a possible trajectory 
control point vector. The number of particles equals 100 for all simulations. For 
larger Z, each particle comprises more dimensions and therefore each particle 
dimension is sampled less densely. 


In order to keep the sampling density in each dimension constant, the number 
of particles needs to be increased exponentially with 7. Since the computational 
effort increases linearly with the number of particles, this quickly becomes 
infeasible especially with computation time restrictions. From these consid- 
erations and the simulation results, it seems beneficial to only use J = lin 
combination with MPF and power penalization. Section 5.3 investigated the 
effect of increasing / in combination with the KF. 


Furthermore, with / = 3 there is less potential for an effect from varying weight 
Rn. Moreover, the approximation for Atpim try = 2 S, Af, =1.5sand/=1 
comprises a larger range of average velocities and energy consumptions than 
the corresponding approximations for At, = 3 s. 
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Temporal safety margin to upper speed limit 


When the temporal safety margin to upper speed limit Ar ;m,rıy is decreased, 
the trajectory follows the speed limit from map data with desired acceleration 
VLim,Map,» more closely which increases the average verlocity. Furthermore, the 
trajectory can follow short lasting inefficient velocity peaks which increases 
energy consumption. Both these effects can be observed for the setting Ar, = 
3 s and J = 3, which achieve comparatively low average velocities. For the 
other settings, however, the average velocity remains almost unchanged but the 
energy consumption decreases. Lowering AtLim,tyy also increases the trajectory 
velocity at local minima of V_im Map,» and can cause the trajectory to exceed 
VLim,Map,y at its local minima. Such a behavior decreases energy consumption 
and increases average velocity at the same time at the expense of driving safety. 
The increase in average velocity is not observed though. 


VLim,Map,» 18 an upper speed limit calculated using parameters that are deter- 
mined in real test drives. Slightly exceeding vLim.Map,» is in most cases not 
noticed but should be avoided because the extent to which the trajectory ex- 
ceeds VLim,Map,» 18 not clearly defined. Figures in Chapter 5 demonstrate that 
trajectories with AfLim,rıy = 1 s stay below VLim,Map,v. Figure 5.4 indicates that 
even with the combination Aflimtry = 1 s, Z = 5 and At, = 5 s the trajectory 


remains below vLim,Map,»- 


For a more complete assessment of the extent and relative frequency of VLim,Map,v 
exceedances, Figure 6.7 shows histograms of the relative difference between 


the actual vehicle velocity vyhe and vLim,Map,» for the AfLim,rjy = 1 s and 


ae 
5000 * 


The histograms reveal that there are exceedances which was to be expected 


Atyimrıy = 2 s with the other parameters Af, = 1.5 s, Z = 1 and Ke = 


because of the unconstrained trajectory optimization problem formulation. The 
majority of these remain below 2 %, which also includes the controller inac- 
curacies. Higher exceedances can be observed but their relative frequencies 
are too low to lead to a safety critical situation during the roughly 1400 s long 
drive. It also needs to be mentioned that map data itself is not perfect and can 
be deprecated or be unavailable, especially during the first seconds after the 
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Figure 6.7: Histograms of the relative frequencies of relative differences between vehicle velocity 
Vvhe! and speed limit from map data with desired acceleration vLim.Map.» for different 
temporal safety margin to upper speed limit At] jm,rıy with RP = sm: Atk = 1.58 
and / = 1. 
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Figure 6.8: Histograms of the relative frequencies of relative differences between vehicle velocity 
Vyncı and speed limit from map data with desired acceleration vLim,Map,v for different 


weight of power error square Rp with Atuimtyy = 15, At, = 1.5 s and 7 = 1. 


driver has taken a different route than expected by the system in real test drives, 
which poses a higher risk in practice. Nevertheless, a slight shift to higher 
velocities is noticeable when Arm rıy is reduced from 2 s to 1 s, which leads 
to a reduction in the energy consumption if the exceedance occurs in curves. 


Figure 6.8 shows histograms of the relative exceedance for two values of the 
weight of power error square Re whereby Atrim try = 15, Af, = 1.5 s and 
I = 1. The results indicate that stronger penalization of the electrical traction 
power mainly leads to an in general slower drive and is unlikely to cause an 
exceedance of vLim.Map,v in curves for the chosen settings. 
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Saving potential and assessment 


O For the parameter combination Atrimtjy = 1 s, At, = 15 s and / = 1, 

the ALC achieves on average the largest absolute energy savings with respect 

to MLC. With additionally Re = m: the ALC requires 3.4 % less energy 

than the MLC at the same average velocity of 59.0 km/h, and for the same 

energy of 16.9 kWh/100km the ALC achieves a 2.6 % higher average velocity. 
1 


. en 1 . . 
Increasing Rp from ¿ppp to zg» reduces the energy consumption by 2.9 % while 


the average velocity decreases by 3.1 %. 


Table 2.2 on page 11 supports assessing the time and energy savings by summa- 
rizing the results for different types of vehicles, roads and ALCs. The results 
in the first three lines were obtained for the WR, however for a conventional 
vehicle. The saving potential of a BEV is much lower than for a conventional 
or hybrid vehicle because of three reasons: 


e Recuperation capability: The table indicates that when investigations 
include a highway section, the saving potential is in general less than 
when only urban and country road sections are considered. This seems 
reasonable as drives on the latter types of roads usually incorporate much 
more acceleration and deceleration than highways and kinetic energy 
cannot be recovered for later use without losses. While conventional 
vehicles need to convert all kinetic energy into heat, hybrid vehicles can 
recuperate at least some of it. BEVs have the highest recuperation power 
limit and therefore can reuse much more energy. 


Efficient power train: The tank to wheel efficiency (c.f. Section 4.4) of 
a BEV is much higher, especially because the electric motor does not 
perform the inefficient conversion from the chemical energy of the fuel 
into kinetic energy as a combustion engine in a conventional vehicle or 
hybrid electric vehicle (HEV) does. An efficiency of up to 95 % for 
electric motors and of up to 35% for petrol engines are stated in [10]. 
As part of the electricity comes from non-renewable energy sources, an 
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inefficient energy conversion process usually also takes place for the oper- 
ation of a BEV but it occurs earlier during the electricity production from 
non-renewable energy sources and is therefore less frequently considered. 
Since the mix of available energy sources varies with the region so does 
the efficiency of fuel and electricity generation. If these processes are 
efficient, the total efficiency also referred to as well to wheel efficiency 
can be twice as high for a BEV than for a conventional vehicle. The gear 
box in a BEV is also more efficient since it usually only has a fixed gear 
ratio and therefore much less moving parts [10]. 


Less degrees of freedom for optimization of driving stategy: Since most 
BEVs have a fixed gear ratio the only degree of freedom is the choice of 
torque, which is however strongly coupled to the trip time and constrained 
quantities such as velocity and acceleration. Opposed to that, the power 
train of a conventional vehicle offers selected gear and clutch state as 
additional degrees of freedom for optimization. In aHEV even the way 
internal combustion engine, electric motor and battery work together can 
be chosen. For that reason the latter vehicle types benefit strongly from a 
look-ahead driving strategy. 


Due to these aspects it is no surprise that the largest time saving of 21.13 % at 
the same energy consumption is stated for a HEV and the largest energy saving 
of 23.8 %, that however is accompanied by a 1.5 % longer trip time, is reported 
for a conventional vehicle. These large savings were achieved in investigations 
including country road and urban sections though. 


For a BEV the results of Table 2.2 state average energy savings of 1.21 % on 
a 1400 km long highway section and of 5.6 % on an only 400 m artificial road 
profile which comes closest to an urban section due to its velocity of roughly 
30 km/h and maximum absolute slopes of 10 %. Compared to these results, 
the savings achieved by the ALC of this work seem reasonable. Adopted from 
[79]. 
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6 Automated energy-efficient longitudinal control 


6.4 Technical contribution 


On the level of driver assistance systems, the contribution of this work is the 
development of an energy-efficient automated longitudinal control for a BEV 
and its parametrization and test in simulations and real test drives. The devel- 
oped ALC system shares some similiarities with the solution presented in [114, 
139, 183]. All systems use predictive route data and measurements of a radar 
sensor. Each system performs an ACC functionality that is extended by the use 
of map data. This enables the system to choose a suitable velocity depending 
on legal speed limits and road geometry. 


The systems can provide the most benefit on winding country roads, on which 
the standard ACC requires frequent manual adaption of the desired velocity. 
The systems can also be used on highways and within cities. However, on 
highways there is a comfort advantage over ACC only when the legal speed 
limit changes and within cities the driver frequently has to interven because the 
systems do no consider right of way or traffic signs. Furthermore, other traffic 
participants are only detected by a radar sensor. 


Power train models enable the systems to determine a more energy efficient 
driving strategy compared to manual driving. 


The differences and unique features of the system described in this work refer to 
the target vehicle and especially the approach for planning the course of velocity. 
The two target vehicles in [139] have an internal combustion engine and the two 
target vehicles in [183] a hybrid power train. Therefore rather complex planning 
approaches that consider gear selection are applied to derive an energy-efficient 
course of velocity. O In contrast a BEV like the target vehicle for this work 
usually has a 1-speed gear box, i.e. a constant gear ratio. The proposed ALC 
includes a novel trajectory planning method that takes advantage of the simple 
BEV power train. 


Compared to the ALC that was developed for the same research vehicle during 
the preceding project, the ALC in this work uses a trajectory with respect to 
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time. Especially at low speeds the time-based approach enables a much more 
dynamic vehicle behavior, which is beneficial for driver acceptance. Further- 
more, the ALC of this work considers the required electrical power explicitly 
and allows to optimize the trajectory not only with respect to travel time and 
comfort but also with respect to driving efficiency. 


In simulative investigations of the energy-saving potential for aBEV the ALC 
requires up to 3.4 % less energy than the MLC for the same average velocity 
and achieves a 2.6 % higher average velocity for the same energy consumption 
on the reference route. 


In a final acceptance test about ten employees of the project partners who all 
work in automotive research and engineering tested the ALC on self chosen 
routes. The ALC yielded throughout good comments, which especially ad- 
dressed the comfortable and smooth driving style. 


The developed system can be used on its own or in combination with other 
assistance systems such as an automated lateral control. Adopted from [79]. 
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7 Conclusions 


7.1 Summary 


This work describes novel methods for the general problem of approximating 
an unbounded number of data points using a B-spline function in the linear 
and nonlinear weighted least squares (WLS) sense. The developed methods 
are based on iterative algorithms for state estimation and their computational 
effort increases linearly with the number of data points. The methods can adjust 
the bounded definition range of a B-spline function during run-time if this is 
required to approximate the latest data point. 


The approximation problem is reformulated as a trajectory optimization prob- 
lem such that the approximation methods compute a velocity trajectory with 
respect to time using data points created from map data. The developed tra- 
jectory optimization methods fall into the category of direct methods and its 
effort increases linearly with the temporal length of the planned trajectory. The 
combination with an adaptive model that describes the power train properties of 
the battery electric vehicle (BEV), allows to plan velocity trajectories whose re- 
sulting energy consumption varies depending on the chosen relative weighting 
of different target criteria. 


The trajectory optimization is extented to an assistance system for automated 
longitudinal control (ALC) that is tested in simulation as well as in real test 
drives. In simulations on a reference route the developed ALC is compared to 
a recorded and re-simulated real drive with manual longitudinal control (MLC) 
and can achieve better ratios between average velocity and energy consumption 
than MLC. 
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7.2 Contribution 


The contribution of this work is threefold and comprises approximation algo- 
rithms, a trajectory planning method and a system for energy-efficient ALC of 
aBEV. 


Iterative algorithms are proposed for the approximation of an unbounded set of 
data as it occurs in online applications. The first algorithm denoted recursive 
B-spline approximation (RBA) and prepublished in [80] solves a linear WLS 
approximation problem iteratively using a Kalman filter (KF). The second al- 
gorithm is termed nonlinear recursive B-spline approximation (NRBA) and 
prepublished in [78]. NRBA is generalization of RBA for nonlinear weighted 
least squares (NWLS) problems and uses a marginalized particle filter (MPF). 
In the MPF a particle filter (PF) deals with the nonlinear subproblem, whereas 
a linear KF solves the linear subproblem optimally. 


RBA and NRBA include a novel shift operation which allows to shift the 
bounded definition range of a B-spline function during run-time such that it 
is always possible to take into account the latest data point. With previous 
recursive algorithms the approximation interval is fixed during run-time and if 
data outside this interval occurs, it cannot be considered. Moreover, the shift 
operation allows to decrease the size of vectors in the filters in order to reduce 
computational effort in both offline and online applications. In offline applica- 
tions all data points are available at once, so their number is bounded and batch 
processing in one step is possible. In contrast, in online applications, additional 
data points are oberved in each time step, therefore a previously calculated 
solution must be updated with each new observation. 


The exponential growth of computational effort with increasing time horizon 
limits the application of direct trajectory optimization approaches to short time 
horizons. O This work presents a direct trajectory optimization method whose 


computational effort only grows linearly with the number of function parame- 
ters or the time horizon. 
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This substantial saving in computational effort is achieved by taking adantage 
of the comparatively simple power train of a BEV with fixed gear ratio and by 
formulating constraints coming from the vehicle dynamics or the environment 
mainly as soft constraints. The resulting nonlinear optimization problem is for- 
mulated as an NWLS problem and solved using NRBA. The trajectory can be 
optimized with respect to travel time, comfort and energy consumption. With- 
out consideration of energy consumption the nonlinear optimization problem 
can be approximated by a quadratic optimization problem that is solved with 
less computational effort using RBA. Adopted from [79]. 


The trajectory is defined by a B-spline function that describes the desired vehi- 
cle velocity with respect to time. The temporal length of the trajectory increases 
with the iterations and the number of function parameters does not need to be 
bounded. If the temporal length of the trajectory must be increased, function 
parameters can be added during the optimization. After each iteration an inter- 
mediate result is available. The optimization can be paused and continued in 
accordance with computation time constraints. 


O Based on the novel trajectory planning approach a system for energy-efficient 
automated longitudinal control of a BEV is developed and tested in simulation 
and real test drives. The system can provide the most benefit on winding country 
roads, but can also be used on highways and within cities. 


In simulations on a chosen reference route the ALC needs up to 3.4 % less 
energy than the MLC for the same average velocity and achieves a 2.6 % higher 
average velocity for the same energy consumption as with MLC. 


For additional subjective impressions, employees in automotive research and 
engineering tested the ALC on other, individually chosen routes. The through- 
out good reviews especially mentioned the comfortable and smooth driving 
style of the ALC. Adopted from [79]. 
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7.3 Suggestions for further research 


The approximation algorithms RBA and NRBA have only been investigated 
for B-spline functions. According to Section 2.3, they can be extended to 
multidimensional control points for approximation of parametric curves and 
surfaces. Furthermore, RBA and NRBA can be transferred to other spline 
representations which offer local control and arise from a linear combination 
of basis functions and control points. 


O Additionally, the modular concept of these algorithms allows to replace the 
included filters. The KF used for the linear approximation is known to be 
optimal, but for the PF and the MPF used for the nonlinear approximation 
various improvements are available in literature. For example, [203] proposes a 
MPF that determines the particles of the PF using particle swarm optimization 
and achieves better results at lower effort compared to both a standard MPF and 
PF. Adopted from [78]. Even more promising for the nonlinear approximation 
seems using a kernel adaptive filter as it does not suffer from the curse of 
dimensionality and requires only solving a convex problem, which is very 
beneficial for real-time applications (c.f. Section 2.4 and Subsection 4.7.2). 


Regarding trajectory planning the computation of the temporal safety margin 
to upper speed limit can be enhanced such that the trajectory follows the upper 
speed limit closer at high speeds while still remaining below local minima of 
the upper speed limit. 


A more comprehensive improvement consists of an adaptive selection of the 
values of temporal distance of neighboring knots, the temporal safety margin to 
upper speed limit and the temporal distance of neighboring data points so that 
the trajectory degrees of freedom per trajectory time unit are not constant but 
depend on the driving situation. For example, on country roads and highways 
with simple shape of the upper speed limit, temporal distance of neighboring 
knots and temporal distance of neighboring data points can be increased. Then, 
for the same computational effort, longer, more farsighted trajectories can be 
planned because the tractory degrees of freedom are temporally less dense. 
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The resulting larger time horizon, that can be optimized simultaneously by the 
proposed local trajectory optimization approach, facilitates achieving higher 
reductions of energy consumption. In contrast, in situations with complicated 
upper speed limit or at low speeds, the values can be reduced to enable a 
more dynamic trajectory. Analogously, the relative weighting of the trajectory 
optimization can be designed adaptive. 


Moreover, additional information can be used for calculating the upper speed 
limit. This can be either data from vehicle-based sensors such as cameras for 
detecting traffic signs and other traffic participants like pedestrians or it can 
be messages of other vehicles or the infrastructure. In particular the driving 
style should be adapted according to the present weather and friction coefficient 
between tire and road [19]. 


The mentioned improvements of the trajectory planning directly translate to 
improvements of the ALC system which can also be extended by other systems 
such as a lateral control system. 


Alternatively, the trajectory optimization approach can be combined with Dy- 
namic Programming (DP), as [187, pp. 1430-1431] suggests and Table 2.5 
illustrates, in order to solve more complex problems, for which the presented 
filter-based approach on its own seems not suitable. 


In such a case, DP can still plan in the spatial domain, which can be even more 
coarsely discretized for less computational effort because a rough long-term 
trajectory is sufficient. From this long-term trajectory the filter-based trajectory 
optimization approach creates a smooth trajectory in time domain with the 
benefits mentioned in Section 5.2 that can represent the vehicle velocity starting 
from zero so that the case of driving off can also be covered by the trajectory 
optimization. 
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